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Prologue 

Recent concern about increased levels of atmospheric CO : and global wanning have led to studies of the 
ocean’s role in the global carbon budget, e.g.. the Joint Global Ocean Flux Study (JGOFS), and increased 
interest in interannual and interdecadal oceanic variability (Ocean Studies Board, [ 1994]; Subcommittee on 
Global Change Research, [1995]). Examination of historical data for these purposes is difficult because of 
the paucity of comprehensive time series data of sufficient accuracy and duration. For meteorological and 
hydrographic data, the ocean weather stations (OWS) and the Comprehensive Ocean- Atmosphere Data Set 
(COADS) [Woodruff et al„ 1987] provide high-quality meteorological observations, and at some OWS 
sites, hydrographic time series are available. However, long time series of biological quantities, e.g., 
profiles of chlorophyll a and column-integrated primary productivity, are relatively scarce, and only 
recently have comprehensive time series been initiated at Hawaii and Bermuda (Karl and Michaels [1996]). 
Alternatively, investigations of the impact of long-term changes in meteorological forcing on ocean physics 
and biology can be pursued using coupled models of physical and biological processes. 

Two contrasting oceanic study areas are addressed in this technical memorandum (TM): the Pacific Warm 
Pool and OWS P (50°N, 145°N). The two areas are different in almost every aspect including seasonality, 
interannual variability, surface flow, and mixed layer dynamics. At OWS P, there is a distinct seasonal cycle 
in all physical and biological parameters, the flow is not well organized and relatively weak, vertical mixing 
is predominantly governed by convective processes, and local physical processes are dominant. In contrast, 
in the Warm Pool region the seasonal cycle is practically absent with interannual signatures (e.g., the El 
Nino Southern Oscillation, ENSO. and remote forcing) being the major contributing factor in the physical 
and biological variability. The flow is very swift and is characterized by the large-scale equatorial current 
system which extends across the entire tropical Pacific. Because of the presence of the Equatorial 
Undercurrent, the vertical changes in velocity are quite intense, and therefore, vertical mixing is 
predominantly controlled by the vertical shear. Also, OWS P is ideal for studying the effects of 
meteorological forcing [Denman and Miyake, 1973; Large et at., 1994] because it is subject to a large 
seasonal temperature (heat flux) cycle and is horizontally uniform with minimal horizontal advection. The 
climatological Ekman upwelling is positive but small because OWS P is located near the climatological 
transition between positive and negative wind stress curl [Reinecker and Ehret, 1988]. Thus, most of the 
vertical fluxes of heat and nutrients are driven by mechanical mixing and buoyancy effects. One important 
conclusion of Martin [1985] was that knowledge of the optical properties that control the vertical 
distribution of light absorption is key to understanding the behavior of the mixed layer and in estimating the 
performance of mixed layer models. For OWS P the optical properties are controlled by biological 
processes and are classified by Mueller and Lang [ 1989]. 

The Pacific Warm Pool is of particular interest because of its role in heat and moisture fluxes between the 
ocean and the atmosphere and, consequently, El Nino and the Southern Oscillation. It was. therefore, the 
site of a special experiment, the Tropical Ocean Global Atmosphere-Coupled Ocean Atmosphere Response 
Experiment, TOGA-COARE, from November 1992 through February 1993 [Webster and Lukas, 1992], 
The longest available time series of observations for equatorial regions are at 165°E. e.g., Delcroix and 
Henin [1991] and McPhaden and Hayes [1991 ]. As part of the Pacific Warm Pool, it is characterized by 
some of the warmest sea surface temperatures (SSTs) and a constant shallow mixed layer (40-50 meters) 
distinctly separate from the thermocline as a result of its low salinity (Lukas and Lindstrom [1991]) which 
presents a barrier to vertical mixing. The data have been extensively analyzed to study the SST variability 
[McPhaden and Hayes, 1991], heat balance [Cronin and McPhaden, 1997], and other ENSO related 
processes (Delcroix et al., [1992]; Delcroix et al. [1993]). Generally, heat budget and circulation studies 
focus on the surface fluxes, but, as noted by Lewis et al. [1990] and Siegel et al. [1995], variations in the 
amount of radiance penetrating through the mixed layer due to phytoplankton absorption can influence the 
heat budget at a significant level. This is particularly true of the Warm Pool where relatively small changes 
in SST and mixed layer heat content can stimulate a substantial response in large scale weather patterns. 
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To address the issues presented above, this TM provides information on the scientific background and 
computational details of a one-dimensional ecological ocean model for equatorial and mid-latitude 
applications. The model can read velocity components (u, v, w) and water temperature (T) from external 
sources, such as field observations or Ocean General Circulation Models (OGCMs), or calculate its own 
ocean temperature profile using one-dimensional mixed layer models (MLMs). There are two choices for 
the MLMs: the Chen MLM [Chen et al., 1994] and the Garwood MLM [Garwood, 1977]. In general, the 
Chen MLM is used in oceanic areas where vertical mixing is predominantly controlled by vertical shear, 
such as the Equatorial Pacific. Conversely, the Garwood MLM was successfully applied to mid-latitude 
areas, such as OWS P (50°N and 165°E), where vertical mixing is predominantly controlled by convective 
processes. 

The following sections provide an overview of the modeling strategy, and a description of the model 
configuration and parameterization to make the user familiar with the model code, pre- and post -processing 
tools, and graphics interfaces that generalize the model applications. 


1. Overview of Modeling Strategy 

1.1 Atmospheric and Marine Optical Models 

As will be described in the following sections, the ecosystem model requires solar irradiance forcing for the 
physical and biological components. Two methods are described here, the Frouin et al. [1989] model and 
the Gregg and Carder [1990] model. The descriptions of these two models, and their respective adaptations 
to the applications at hand, were reproduced from McClain et al. [1996] and McClain et al. [1998, 
submitted], respectively. 

1.1.1 Non-Spectral, Clear Sky Scheme 

Frouin et al. [1989] estimated the downwelling clear sky solar irradiance £ deilr for two spectral domains, 
350-700 nm (photosynthetically available radiation, PAR) and 250-4000 nm (total spectrum). PAR, rather 
than spectral irradiance, can be used to drive phytoplankton growth in low biomass (Morel and Prieur, 
[1977]) waters where only one phytoplankton species of constant specific absorption is assumed. In fact, 
simulations using a full spectral version of the coupled ecosystem model produced nearly identical monthly 
climatological chlorophyll a and temperature profiles to the OWS P application study, but at significant 
computational expense. 

The near-infrared irradiance is calculated using the approximation 


£c , ear (700 “ 4000) = £ clear (250 - 4000) - £ clear (350 - 700). 


The error introduced by including the irradiance between 250 and 350 nm is small. 

The model input includes ozone concentration, aerosol type (marine and continental), and visibility. A 
constant ozone concentration of 340 Dobson units is used, and the aerosol type is assumed to be marine. 
Visibility (V0 is computed at each time step from the relative humidity (for the range of 65% to 99%) using 
the equation (Neuberger, [1951]): 


V = 231 -2.3 RH , 


(2) 


where RH is the relative humidity (in percent). For RH < 65%, v = 80 km, RH is determined using 
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RH = 


100 e^T) 

**<T) 


(3) 


where e u is the vapor pressure at the observed temperature and t\ ir is the vapor pressure at the observed dew 
point temperature (Smith. [ 1966]). Vapor pressure is related to temperature by 

e(X) = a i ,e l “ 1+T ', (4) 

where a„ = 6.1078. a, = 17.269, a : = 237.29°, and 7 is the temperature (degrees Celsius) of either air or 
water [Tetens, 1930], 

The irradiance that penetrates the ocean surface, E{ O'), which is a function of the clear sky irradiance 
attenuated by cloud effects (see the calculation of £ surf in Section 1.1.3), is computed as 

£(0-) = (l-//)£ surf . (5) 


where |i is the surface reflectivity. The surface reflectivity consists of three components, 

M-Mj+Mf+Mss- (6) 

where |i (/ is the direct specular reflectance, is the reflectance from foam, and \i ss is the subsurface 
reflectance due to molecular and particulate backscatter. The direct reflectance is taken to be the Fresnel 
reflectance for a flat surface unless the solar zenith angle is greater than 40° and the wind speed, VV, is 
greater than 2ms 1 . If these two conditions are satisfied [Austin, 1974], then 

fl t , = 0.0253/ ,e _40 \ (7) 


where 0 O is the solar zenith angle and 


b = -0.000714^+0.0618. (8) 

The foam reflectance is computed as the product of the effective reflectance of foam (0.22 [Koepke, 1984]) 
and the wind-speed-dependent fraction of the surface covered by foam [Monahan and CTMuircheartaigh. 
1980] or 


ju f =6A9x\0~ 1 W X52 . <9) 

For 0.01< chlorophyll a < 0.5 mg m \ p„ is approximately 0.035 [Morel. 1988]. 

After obtaining £(0), the downwelling irradiance at depth r (the vertical coordinate decreases in the 
downward direction) is estimated by 

£(r) = E(z~A:)e ~ K[z) ^ , (10) 


where K is the downwelling diffuse attenuation coefficient. The K associated with each wavelength domain, 
PAR and the near-infrared, is determined independently. £(PAR) is computed as a function of the chi a 
concentration at each depth (Morel [ 1988]; see also Wroblewski and Richman 1 1987] and Frost 1 1993]): 

K( r, PAR ) = K K (PAR) + 0.05 1 8chl(z) ,,42x , (11) 

where K w (PAR) is the effective clear water diffuse attenuation coefficient (0.027 m l ) and chl(r) is the 
chlorophyll a concentration at depth r. It should be noted that the monthly mean profiles of PAR obtained 
from this formulation were not substantially different from those obtained when integrating the spectral 
version of the model from 400 to 700 nm. For the near-infrared, K (NIR) is the effective value of (NIR) 
over the wavelength band of 700-4000 nm and is assumed to be a constant =2.5 m 1 (Baker and Smith, 
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[1982]). The two wavelength bands are handled separately in the heat absorption calculation. The near- 
infrared radiation is essentially captured within the first meter of the water column. 

The rate of heating at a given depth because of the absorption of downwelling irradiance is given by 

dT{z) = {[E(z-Az,PAR)-£UPAR)] + [E(r-Az,NIR)-E(z,NIR)]}Az 

* “ P* c r» “ 2I 

where p M and c pw are the density and heat capacity, respectively, of water. 


1.1.2 Clear Sky Spectral Atmospheric and Marine Optical Schemes 

The model of Gregg and Carder [1990] can be used to compute the clear sky irradiance in the visible from 
280-700 nm. The published version of the model only includes 350-700 nm (PAR), but is easily extended 
to include shorter wavelengths given the solar constants and ozone absorption coefficients for those 
wavelengths. Gregg and Carder’s algorithm requires estimates of surface wind speed (from the European 
Center for Medium Range Weather Forecasting [ECMWF] monthly data), relative humidity (from Bishop 
and Rossow [1991] ), precipitable water (from the International Satellite Cloud Cover Project [ISCCP]), 
visibility (from Bishop and Rossow [ 1991] ) and ozone concentrations (from ISCCP). 


1.1.3 Cloud Cover Corrections 

Three different approaches for incorporating cloud cover effects on the radiative transfer models were used: 
two of these approaches were specifically designed for OWS P and the Warm Pool applications, whereas 
the third approach, developed by Tai and McClain [1996], is only possible when broadband, clear sky 
radiance observations are available for the region of interest and time period of simulation. The algorithm 
of Frouin et al. [1989] is used to calculate the total and PAR radiances for clear skies if observations are 
not available. A four-year (1987-1990) comparison was conducted between the total clear sky radiance 
provided by the Frouin et al. [1989] method and the equivalent clear sky radiance provided by Bishop and 
Rossow [1991]. The 4-year mean values of clear sky radiance derived from these two sources agree within 
3% (295 W/m“ from Frouin et al. [1989] and 303 W/m 2 from Bishop and Rossow[1991]). 

The following is a description of the three different cloud cover corrections. 

OWS P 

This approach consists of an exponential fit to the fractional cloud cover. Cloud cover data for the period 
of January 1951 through December 1980 were obtained from the OWS P surface observation time series. 
These data have a resolution of 3 hours and vary from 6 oktas in autumn to 7.75 oktas in summer. By 
comparison, the COADS monthly mean cloud cover climatology has a similar pattern, with slightly higher 
values in autumn. 

Clear sky irradiance is modified to account for the observed cloud cover by applying a power law correction 
tuned to yield the observed climatological monthly mean surface irradiances [Dobson and Smith, 1988], 

fur, =f,ear0 -0.53c 0 ' 5 ), (13) 

where £ slirf is the downwelling surface irradiance and c is the fractional cloud cover. 


Warm Pool 

Correcting the clear sky values for cloud cover is more difficult. For example, in the Warm Pool area, total 
surface irradiance observations from the TAO array are limited to August 1992 through April 1994. Total 
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monthly mean surface irradiance and cloud fraction estimates from the ISCCP [Schiffer and Rossow, 1985] 
are only available for July 1983 through June 1991. ISCCP surface PAR estimates [Bishop and Rossow, 
1991 1 are not available for the time period of interest ( 1983-1992). Therefore, an algorithm for estimating 
the visible surface spectral irradiance from the theoretical clear sky irradiances and monthly mean ISCCP 
total irradiances and cloud cover was developed. 

The approach emulated the ISCCP processing by using a proxy formulation which incorporates coefficients 
derived from statistical comparisons with ISCCP data. ISCCP total (tot) surface irradiance (250-4000 nm) 
and PAR are approximated using 


£ (tot) = E c f (tot)Cj , 

(14) 

£ (PAR) = E c f (PAR XT, , 

(15) 


where (tot/P AR) is the equivalent ISCCP irradiance, E F (tot/PAR) is the clear sky irradiance computed 

using the Frouin et al. [1989] algorithm, and C, is the equivalent ISCCP cloud cover correction. The 
Frouin et al. [1989] algorithm requires water vapor (precipitable water), visibility, and ozone input. A 
marine haze model is assumed to be independent of wavelength, at least to the first order. 

To estimate C, using ISCCP cloud fraction data, the cloud correction scheme of Dobson and Smith [1988] 
is adopted, i.e.. 


Eds ( tot ) = E o (tot) [cos 0 O (A, + B, cos 6>J] , (16) 

where E DS (tot) is the Dobson and Smith [1988] total irradiance, E„ is the extraterrestrial irradiance (1358.2 
W m 2 ), 0 O is the solar zenith angle, and A, and B i are empirical constants for each cloud okta (i.e., 0-8, okta 
= 0 will be designated as “c” for clear sky). Dividing (16) by (14) yields 


Bps = E PS PS 

£,(tot) E ( r {ioi)C f 


where 


£’p 5 (tot) = £„(tot)[cos0„(A c +B C cos0J] 


(18) 


and 


A t + B i cos0 o 
A c + B c cos0 o 


(19) 


The values for the clear sky coefficients A c and B, are 0.40 and 0.386, respectively. The values for A , and B, 
vary according to the oktas. 

Comparison of monthly mean values of E DS computed using monthly mean cloud cover with E, yielded a 
mean value of 1 . 1 1 for the left side of (17). Similarly, an annual cycle of E L F compared with E DS values 
produced a fairly constant ratio of 0.95 for the first term on the right side of (17). It is assumed that these 
ratios are independent of wavelength and can be applied in the visible range. 
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Frouin et al. [1989) provided algorithms for both PAR and total surface irradiance while Gregg and Carder 
(1990) only addressed spectral PAR. Finally, the spectra for PAR can be expressed as 

E,U) = 0.85/( par e£ c (A)C ds (20) 

where R PAR = E F (PAR)/ E 6r (PAR) at each time step. No smoothing of the monthly mean cloud fraction 
was included in the model surface irradiance computation. The 4-year mean (1987-1990) cloudy sky 
radiance derived from (14) is 243 W/m 2 . This value is within less than 1% of the equivalent 4-year mean 
cloudy sky radiance (241 W/m 2 ) obtained from the Bishop and Rossow [1991] data set. 

Cloud Attenuation Scheme usins Observed Irradiances 

When broad band in situ irradiance data are available, a more accurate scheme for cloud attenuation can be 
used. Leonard and McClain [1998, submitted] applied the scheme of Tai and McClain [1996], which is 
based on a combination of broad band irradiances derived from both irradiance models and field 
observations. The scheme uses a ratio between the broad band irradiance from the Tropical Atmosphere 
Ocean (TAO) buoy at 140°W on the equator (£obs)' and the irradiance (E F c ) calculated using the Frouin et 
al. [1989] broad band irradiance model. The modeled PAR was then multiplied by that ratio to obtain 
cloud-corrected PAR. Specifically, 


E 


PAR 



/ ^OBs(tOt) 

E c f (tot) 


(21) 


The modeled, cloud-corrected surface PAR compared well with in situ PAR data from the Joint Global 
Ocean Flux Study (JGOFS) equatorial Pacific cruises in 1992, validating the assumption that the 
attenuation of light by clouds is not strongly wavelength dependent. 

1.2 The Garwood MLM 


As mentioned in the Prologue of this TM, the Garwood [1977] mixed layer model was used to simulate 
mixed layer dynamics for the OWS P study because it performs best in oceanic areas where mixing is 
predominantly controlled by vertical convection, e.g., mid latitudes. 

Garwood [1977] developed a one-dimensional bulk model of the mixed layer of the upper ocean. The 
model includes an entrainment scheme dependent on the relative distribution of turbulent energy between 
horizontal and vertical components as a mechanism for both entrainment and layer retreat. The model has 
two key dynamical properties for the generation and dissipation of turbulent energy: 

• The fraction of wind-generated turbulent kinetic energy partitioned to potential energy increase by 

means of mixed layer deepening is dependent upon layer stability, H* = /i/L, as measured by the ratio 

of mixed layer depth h to Obukhov length L. This results in a modulation of the mean entrainment rate 
by the diurnal heating and cooling cycle. 

• Viscous dissipation is enhanced for increased values of the inverted Rossby number, R 0 1 =/?/*/ w*, 

where f is the Coriolis parameter and u * the friction velocity for the water. This enables a cyclical 

steady state to occur over an annual period by limiting maximum layer depth. 


The main model assumption is that the turbulence of the overlying mixed layer provides the energy needed 
to de-stabilize and erode the underlying stable water mass. The turbulent kinetic energy budget is, 
therefore, the basis for the entrainment hypothesis: 
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1 d(u~ + V + M’ j p 

2 dt "~ 


dU 

UW 1- I’M’— — 

ct dz 


1 — d 

( 

+ bw 

w 

J A 

L ^ 


U~ + A" + M’~ p 
+ — 

2 P„ 


— ~ 0, (22) 


where € is the viscous dissipation, p is pressure, and the ensemble mean and fluctuating components are 
denoted by the upper and lower cases, respectively. If the terms in (21) are known throughout the boundary 
layer, then the evolution of the potential energy and density profile may be evaluated using the budget for 
mechanical energy. The buoyancy equation is generated from the heat and salt equations together with an 
equation of state 

p = p^-a(e-0 c )+p(S-S,)\, (23) 

and the definition for buoyancy 


b = g(p„-p)/P<, • (24) 

In (22) and (23) 0 is the potential temperature, S salinity and p the density, while a and p are the expansion 
coefficients for heat and salt, respectively, and g is the gravity. 

The bulk equations to be solved are the vertically integrated (within the mixed layer depth h) equations for 
buoyancy and momentum 


+ (25) 

dt dt p o C p v ' 

/j + A = fh(v ) - wvv(O) , and (26) 

ot ot 

h^^ + AV^-A = -fh{U}-™(0), (27) 

ot ot 

where Q a is the surface heat flux, C p is the specific heat of sea water at constant pressure, and A is the 
Heaviside step function. B. U, and V refer to the time averaged values of buoyancy and velocity 
components. The coordinate system has x positive to the east, v positive to the north and r positive upward. 
Three assumptions were made in this integration: 

(i) Vertical fluxes are negligible below the mixed layer, therefore, 

hw(-h - 5) - uw(-h -S)= vw(-h -S) = 0 . 


(ii) The mixed layer is sufficiently homogeneous so that 

A B~(B)- B(-h - 8 ) 
A U ~ (il)-U(-h - 8) 
A V ~(V)-V(~h -8) 


(iii) Horizontal homogeneity is assumed for all mean variables. 
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1.3 The Chen et al. [1994] MLM 


The Chen et al. [1994] MLM is based jointly on the Kraus-Tumer type mixed layer model and Price et al. 
[1986] dynamical instability model. The scheme is computationally efficient and is capable of simulating 
the three major mechanisms of vertical turbulent mixing in the upper ocean, e.g., wind stirring, shear 
instability, and convective overturning. 

The model ocean is divided into a fixed number of layers. The first layer represents the mixed layer, and 
the layers below the mixed layer are chosen according to a a coordinate. In this fashion, the thickness of all 
layers is allowed to change, but the ratio of each a layer to the total water column below the mixed layer is 
fixed to its prescribed value. Once the mixed layer depth is predicted, the thickness of the remaining layers 
are calculated diagnostically. The real advantage of such vertical coordinate system is its computational 
efficiency. Finer vertical resolution is assigned where it is needed, i.e., below the mixed layer [Gent and 
Cane, 1989]. This overcomes one of the disadvantages of the level models which need high resolution over 
the whole range of the mixed layer variation. 

Since the layer interfaces are not material surfaces, there are fluxes of mass, buoyancy, and momentum 
through these interfaces. Of particular importance is the interface between the mixed layer and the layer 
right below, because the fluxes through it are directly related to the atmospheric forcing at the ocean 
surface. 


The one-dimensional model equations are 


dt 


(h k u k )+(wu \_ [/2 -(wu ) i+1/2 + Jkxh t u t = rj.,,, - r t+l/2 , 


(28) 


^(*A)+( W *)m/2 - (W>), 


*+ 1/2 


~ ^k-\n 


*+ 1/2 > 


( 29 ) 


dh k 

+ W *- l /2 ~ ^*+ 1/2 = 0 ( 30 ) 

k = 1,2, ,nz 

where u is the velocity vector, /is the Coriolis parameter, k is the unit vertical vector, t is the wind stress, 
and B is the buoyancy flux. For simplicity, buoyancy b is assumed to be a linear function of temperature T . 
Thus (28) can be rewritten as 

) + _ ( W ^)*+!/2 = Qk~ 1/2 ~Qk+ 1/2 ( 31 ) 

where Q is the heat flux. A similar equation for salinity, S , can be included and b is determined from T and 
S using a simplified equation of state [Mellor, 1991], Generally, t and Q take the form 


( 32 ) 


^A+l/2 


JU1/2 ( U A U * + l ) 


^ _ r^K + 1/2 (7* ^ / 

AT + 1/2 ** T ; -h/ 4+!/2’ 

/? 4 


where A'v/ and K r are vertical eddy viscosity and diffusivity, respectively, and I-rl n e K ~ is the penetrating 
solar radiation, where r is the fraction of the total irradiance, /„, that is in the ultraviolet-visible portion of 
the spectrum, and K is the average diffuse attenuation coefficient for the ultraviolet-visible spectrum over 
the mixed layer. K and r are assumed to be constant by Chen et al. However, both are functions of time 
and K is a function of depth as well. The ecosystem model computes both r and K at each time step and the 
values can be provided to the mixed layer model in the coupled mode. At the top of the model, T 1/2 and Q 1/2 
are specified as surface wind stress and heat flux; at the bottom of the model, t nz+ j /2 and Q^+uz are set to 
zero. The coefficients K M and K T are set to small constant values of background mixing {Km -1 and K r = 
0.1 cm 2 s' 1 ). The principal turbulent mixing is caused by mixed layer entrainment, convection, and shear 
instability. The role of entrainment or detrainment is parameterized in (27) and (30) by the terms 
containing w. The convection and shear-produced mixing, however, are not explicitly included in (27) and 
(30); they are instantaneous convective adjustments. 

1.3.1 Heat Budget Adjustment for the Chen et al. MLM 


While simulating temperature profiles with the Chen et al. [1994] mixed layer model coupled with an 
ecosystem model (discussed below) in the Warm Pool, a slow warming trend was noticeable at all layers (0 
to 250 meters). This warming trend was monotonic at about 6.3 x 10 s °C s' 1 , or 8°C in 4 years (the length 
of the simulation). Figure 1 shows the time series of T(z); notice the downward trend of the isotherms 
towards the end of the simulation. Figure 2 shows the equivalent time series of T(z) obtained from TOGA- 
TAO data at exactly the same location and time period. A comparison of Figures 1 and 2 reveals that the 
warming trend is obviously anomalous. 


The authors attributed this anomalous wanning trend to the limiting dynamics provided by the one- 
dimensional mixed layer model. The heat balance in the equatorial region cannot be accurately modeled 
without considering horizontal advection and diffusion terms. The one dimensional MLM does not 
exchange heat through the bottom of the vertical domain (250 meters). Therefore, the excess heat that 
would be exchanged laterally via the advection and diffusion terms in the three-dimensional 
thermodynamics, is constantly being accumulated in the one-dimensional vertical layers of the MLM. 


To compensate for the anomalous warming, a strategy was developed to allow the excess heat to be flushed 
out from all layers of the MLM. This is an ad hoc but necessary strategy to ensure that an accurate T{z) is 
provided to the ecosystem model. This strategy consists of an explicit finite difference for the heat equation 
that is applied at the end of each time step to T{z) calculated implicitly by the original MLM 
thermodynamics. Specifically, 


3T_ 

dt 


- -w 


dr 

— + 

* 



V 


\ 

/ 


is subject to the following boundary conditions: 

(dr\ 


\ ^ j --= 0 


= o, 

T(-h)=T hil , ll>m = \rc. 


(33) 


( 34 ) 
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Figure 1. Four-year (1987-1990) time series of monthly-averaged MLM temperature at 165°E on the 
equator. The heat flux correction was not applied. 
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Figure 2. Four-year (1987-1990) time series of monthly-averaged TOGA-TAO temperature at 165°E on 
the equator. 
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The finite difference form of (32), which is up-wind in time and centered -difference in depth, can be 
written as 


/j+i 


r n . . 

- W 


r r^ ft nr n 

l k ~ 1 U\ 

A 2 


A t + K. 


*-p n 
1 *-l 


A - 2 


.? T" 

— -At 


(35) 


where n is the time level. At is the time step, k is the index for the vertical level, and Ar is the layer 
thickness. Numerous sensitivity runs were conducted with different values of K v to find the required value 
that compensates for the anomalous warming. The best value for K v was 3.3 x 10' 4 m 2 s' 1 (3.3 cm 2 s' 1 ). 
Figure 3 shows a time series of T(z) after the heat flux correction has been applied. Note that the isotherms 
are now in much closer agreement with the TOGA-TAO time series shown in Figure 2. Also, in Figure 4, 
a comparison is shown between the time series of T at 250 meters from the MLM and the time series of T at 
250 meters from TOGA-TAO, before and after the heat flux correction (top tier and lower tier, 
respectively). Note that the warming trend has been completely removed from the record once the heat flux 
adjustment has been applied. 

To verify how realistic this heat flux would be in terms of an equivalent horizontal flux, the required zonal 
temperature gradient for heat advection, udTfdx , is calculated and the result with values of zonal 
temperature gradients derived from TOGA-TAO data are compared. At 250 meters, the zonal component 
of the velocity, «, averages 20 cm s 1 (from TOGA-TAO data). With this value of w, the required 3775* 
that matches the adjustment value of 6.3 x 10' 8 °C s' 1 is 3.2 x 10' 7 °C m'\ or 0.032 °C/100 km. According 
to Wang and McPhaden [1998], zonal temperature gradients of this magnitude are easily found in the 
Warm Pool area and, therefore, the hypothesis of heat leakage via horizontal advection is a feasible 
thermodynamic mechanism in the heat budget of the western equatorial area. 

1.4 The Ocean General Circulation Model (OGCM) 

The OGCM is the reduced gravity, primitive equation, a coordinate model of Gent and Cane [1989] with 
an embedded hybrid mixing scheme of Chen et al. [1994]. Surface heat fluxes are computed by coupling 
the OGCM to an advective atmospheric mixed layer model (Murtugudde et al., [1996]). Complete 
hydrology has been added to the model with freshwater forcing treated as a natural boundary condition 
(Huang, [1993]). The UNESCO equation of state is used for computing buoyancy from salinity and 
temperature. The improvements in tropical SST simulations and upper ocean hydrology was reported in 
Murtugudde et al. [1995] and Murtugudde and Busalacchi [1997], respectively. The vertical structure of 
the model ocean consists of a mixed layer and a specified number of layers below according to a a 
coordinate (total number of layers is 20 for the Warm Pool simulation). The horizontal grid is stretched 
down to 1/3° resolution near the equator and at the eastern and western boundaries. The model domain 
spans the Pacific zonally with meridional boundaries at ±30° latitude. The mixed layer depth and the 
thickness of the last a layer are computed prognosticate and the remaining layers are computed 
diagnostically such that the ratio of each o layer to the total depth below the mixed layer is held to its 
prescribed value. The model has a 1-hour time step and 10-day average velocity, temperature, and mixed 
layer values are output tor use in ecosystem models with one-way coupling. The model is spun up with 
climatological winds for 10 years. The interannual simulation for 1983-1992 is initialized with the 
climatological run and forced with the monthly mean Florida State University (FSU) winds. The 
precipitation data of Xie and Arkin [1996] was used for freshwater forcing. 

1.5 One-Way Coupling with Ocean Global Circulation Model 

The ecological model can be driven by velocity («,v, w) and water temperature fields (T) originating from 
the OGCM output. In this manner, the OGCM is configured and runs independently from the ecological 
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Temperature in °C (Kv = 0.000333 m 2 /s) 



Time in Months (01/87 to 12/90) 


Figure 3. Four-year (1987-1990) time series of monthly-averaged MLM temperature at 165°E on the 
equator. The heat flux correction was applied. 
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Temperature Time Series without Flux Correction 



Temperature Time Series with Flux Correction (Kv = 0.000333 m 2 /s) 



Figure 4 . Four-year (1987-1990) daily time series of TOGA-TAO and MLM temperature at 250 m. The 
upper tear shows the comparison between the MLM temperature without heat flux correction and the 
TOGA-TAO temperature. The lower tear compares the MLM temperature using the heat flux correction 
and the same TOGA-TAO temperature time series. Note the elimination of the warming trend in this lower 
tear. 
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model; the «, v, w, and T fields are stored at preset time intervals (in this case, 10 day averages) for 
subsequent use as input to the ecological model The ecosystem model is executed at one hour intervals 
with OCGM input derived from a linear interpolation of the 10-day average values. Other combinations for 
the driving physical fields have also been tried. An example is the usage of the vertical velocity component 
from the OGCM and w, v\ and T fields from the TOGA moorings. Figure 5 shows the processing flow 
when the ecosystem model is linked to either a mixed layer model (coupled mode) or is forced by 
prescribed input (uncoupled or 1-way forcing mode). 


1.6 Ecosystem Model 

The ecosystem model (Figure 5) is essentially the same as that used in McClain et al. [1996] with the 
following modifications: 

1. The light model used in the Warm Pool and OWS P applications includes spectral absorption from 280 
to 700 nm. The spectral light model requires absorption spectra for water and chlorophyll a . The 
seawater absorption spectrum of Baker and Smith [1982| was used. The use of a spectral model 
eliminated the need to estimate light attenuation for PAR, i.e., XXPAR), using the Morel [1988] 
empirical function based on chlorophyll-a concentration. Also, the growth rate and C:Chl (carbon to 
chlorophyll a) ratio computations require integration over PAR in units of pmol quanta. The 
transformation from W m 2 nm 1 to pmol quanta m 2 s' 1 nm 1 was applied to each wavelength according 
to the expression (Zeebe et al., [1996]) 


E Mm '“(A) = 


£ w / UP 6 
N a hc 


(36) 


where /V A is the Avogadro number (6.022. 10 21 ), h is Planck’s constant (6.626. 10 34 J s), c is the speed ot 
light (2.998.10 s m s' 1 ), and X is in meters. 

2. In the Warm Pool application, the surface-wind-stress-dependent algorithm for the eddy diffusion, K v , 
has been replaced with the method of Pacanowski and Philander [1981], which uses profiles of 
temperature and horizontal currents ( u and v) obtained front the OGCM. 

3. In the Warm Pool application the vertical velocity, n\ is obtained from the OGCM rather than being 
derived from the curl of the local wind stress. 

4. In the Warm Pool application the bottom boundary (250 m) nitrate concentration is allowed to change as 
a function of temperature, T (° C). A variable bottom boundary condition was needed because, at the 
equator, planetary waves introduce significant vertical excursions of the isopycnal surfaces below 

250 m. Using climatological nutrient data at 165°E from Conkright et al. [19941 for standard depths 
between 150 and 350 m. 


NO, =39.64-1.4677' 


(37) 


with an r 2 = 0.92. 


5. For simulations in the Warm Pool, the exponential temperature-dependent zooplankton respiration rate 
(;-.) algorithm was replaced with a linear relationship based on data from Ikeda 1 1985], The linear 
respiration relationship. 


r. = 0.0722 + 0.0067(7 - 20), 


(38) 
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Figure 5. Flow chart ot ecosystem model showing the processing flow and the linkage between the four 
component biological model with the physical fields. The chart shows two options for the use of the 
physical fields; one originating from the one -dimensional MLM and another originating from the OGCM. 
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produces rates of less than 0.15/day. 

6. In the Warm Pool application, the half saturation constant for NH 4 , was increased from 0. 1 to 0.5. 

Also, the NH 4 inhibition coefficient, pk. was increased from 3 to 5. These changes were introduced to 
help limit growth rates to a range consistent with observed values since maximum growth rates predicted 
by the Eppley [1972] expression are ~ 5 d' 1 for the surface temperatures in the Warm Pool. 

7. For the Warm Pool simulations, phytoplankton mortality, nu was increased from 0 to O.id 1 to obtain 
phytoplankton concentrations more consistent with observations. Also, zooplankton mortality was 
increased from 0.04 to 0.15 d 1 to obtain a pronounced subsurface chlorophyll maximum. 

8. The phytoplankton sinking rate was set to zero to reflect the tact that cell sizes are relatively small in the 
Warm Pool [Blanchot et al., 1992]. 

9. The C:Chl a ratio in the Warm Pool application was limited to the range 20-300 to be consistent with 
the observed range of values [Furuya, 1990; Chavez et al., 1991 ] for the equatorial Pacific. 

10. In both Warm Pool and OWS P applications, fecal pellet remineralization (conversion to NH 4 ) was 
added as a constant fraction (c^) of the fecal pellet production. About 70% of the pellets are recycled 
in the photic layer of the HNLC (high nutrient-low chlorophyll; Dam et al., 1995) and 84% [Le Borgne 
and Rodier, 1996] in the oligotrophic regimes of the equatorial Pacific. A value of = 0.8 was used. 

1 1. For both Warm Pool and OWS P simulations, an ammonium nitrification algorithm based on Olson 
[1981 J was incorporated. The strategy is to model irradiation dosage between 300-470 nm to inhibit 
bacterial conversion of NH 4 to NCE near the surface. The addition of this process was prompted by the 
buildup of NH 4 at depths below the observed subsurface NH 4 maxima between 60-120 m (James 
Murray, personal communication). The rate of conversion of ammonium to nitrate, A n (nmol/day) is 


*■(-’)= aim 


p( 

£>(--)- K.CO 


(39) 


where 


1 7 470 

D(z) = - J J E s (z,A,t).as(A).dM (40) 

* 0 3(H) 

T is 24 hours, D min = 0.0095 W/m 2 , 4" max = 2 nmol/day. K r > = 0.036 W/nr, and E s has units of W/m 2 /nm. K D 
is the half saturation dosage or the dosage at which inhibition is half the maximum rate, and D min is the 
minimum inhibition dosage. The absorption coefficient for NH4, as(X), was digitized from Figure 3 of 
Olson [1981] which provides the action spectra for photoinhibition of Nitrobacter winogradslcvi and 
Nitrosomonas europaea ; the photoinhibition spectrum for the latter species was used. 

The following system of coupled differential equations simulates the dynamics of phytoplankton nitrogen 
(P), micro-zooplankton nitrogen (Z), ammonium (NH 4 ), and nitrate (NO3) stocks within the upper ocean: 


dP dP dSP df., dP\ 


GP - mP — r p P — IZ , 


(41) 
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and 


dZ dZ 




- X)1Z - gZ - r z Z , 


(42) 


r)NH 4 

a 


4- Vi ' 


* 


*1 " * J~ 

( a pm + r P ~^iG)p + (a : g + r.+ c p( ,M )z - A" , (43) 


dNO, 

a 


•+ w. 


dNC > 3 
dz 



dm: 

dz j 


— — Tt^jP + A n . 


(44) 


Table 1 provides the variable definitions and the values of the free parameters used in deriving various 


terms in these equations. The rate of phytoplankton growth (G), the herbivore grazing term (/), and the 
partitioning of NH 4 and N0 4 uptake (7t) are calculated using the following equations: 

G = PG m , 

(45) 

G m =G/*' T , 

(46) 

I=RJ \-e^ p ]. 

(47) 

x = NH - 
' NH 4lim + N0 31im ’ 

(48) 

x _ N0 3.im 

2 NH 41m) +N0 31im ’ 

(49) 


where, in (44) and (45), G m (d ') is the first-order exponential function of temperature (Eppley, [1972]), G„ 
(d ') is the specific growth rate at 0°C, k gl> (°C ' ) is a rate constant which determines the sensitivity of G„, to 
changes in T. (3 is the resource limitation term calculated as the smaller of the computed nutrient and light 
limitation terms (N nm and L Um , respectively). In (46), R„, is the maximum ingestion rate of zooplankton and 
A determines the sensitivity of grazing to phytoplankton stock. Ammonium (NH 4 | im ) and nitrate (N0 3hm ) 
limitation are calculated using the standard hyperbolic formulation (Monod [1942]). More specifically, in 
Equations (47) and (48), the ammonium limitation term is calculated using the equation 


NH 


4 Jim 


nh 4 

*nh 4 +NH 4 


(50) 


where is the NH 4 concentration where growth is half maximal. The nitrate limitation term is 

calculated in a similar manner but with the addition of a term to correct for inhibition of N0 3 uptake due to 
the presence of NH 4 (Wroblewski [1977]) 
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NO 


3 lim 


NO, 


,-/>* nh 4 


^no,+NO, 


(51) 


where A' NO , is the NO, concentration where growth is half-maximal and pk detennines the sensitivity ot NO, 

uptake to the presence of NH 4 . The nutrient limitation term. N, im . is defined as the maximum value between 
NH 4 |j m and NO.,| im . The light limitation term, L lim , is a function of the ambient irradiance and the 
photoadaptation parameter (McClain et al. [1996]). 

The coordinate system origin is at the surface and negative downward. To prevent numerical instability 
when solving the finite difference form of the continuum equations above, the Crank-Nicholson scheme 
(Press et al., [ 1988]) was applied. This involves the simultaneous solution of these linear equations for each 
time step and each depth using the method of Gaussian elimination for tridiagonal systems. The Neumann 
boundary condition, dX/dz = 0, is applied at both the surface and lower (350 m for OWS P and 250 m for 
the Warm Pool) model domain boundaries for P and Z. Initial profiles of temperature and NO, are usually 
obtained from annual climatologies. Depth-independent initial concentrations of P, Z, and NH 4 are given in 
Table 1. For the OWS P study, NH 4 and NO, were given a fixed value equal to the initial condition at the 
lower boundary and the Neumann condition was applied at the surface. 


2. Model Applications 

2.1 Ocean Weather Station P, 1951-1980 

Physical and biological processes in the mixed layer at OWS P (50°N, 145°W) over a 30-year period 
(1951-1980) were investigated using observations and model simulations (McClain et al., [1996]). The 
observations include 30 years of surface meteorological and sea surface temperature data collected at OWS 
P and Ekman upwelling velocities derived from COADS, 14 years (1953-1966) ot daily temperature 
profiles, nearly 150 chlorophyll a profiles spanning all months of the year, monthly climatological solar 
radiances, and 0- to 50-m integrated nitrate concentrations. The simulations incorporated models for the 
estimation of surface solar downwelling irradiance, surface heat fluxes, subsurface diffuse attenuation, 
mixed layer dynamics, and biological processes. The time-dependent model input was the surface 
observations of cloud cover, air temperature, dew point temperature, and wind speed. The atmospheric 
irradiance, marine diffuse attenuation, and mixed layer models were adapted from existing models 
developed by others. 

The mixed layer depth (MLD) is calculated using an updated version of the Garwood [1977] one- 
dimensional dynamical MLM. as described in Section 1.1. This model predicts the rate ot deepening or 
retreat of the mixed layer as a result of entrainment or detrainment processes at its lower boundary. 

The biological model, developed by McClain et al. [1996], has four components (nitrate, ammonium, 
phytoplankton nitrogen, and zooplankton nitrogen) and computes a variety of additional quantities 
including chlorophyll a concentration and gross and new production. The non-spectral model produced 
seasonal climatological chlorophyll a protiles within the observed standard deviations at almost all depths, 
but with values markedly less than the mean observed values in some seasons, e.g., summer. Also, the non- 
spectral model did not include nitrification of NH 4 to N0 3 resulting in the accumulation of NH 4 at depth. 
Climatological monthly profiles of temperature were within one standard deviation ot the observed values 
at nearly all depths. Also, the climatological annual cycles of solar irradiance and 0- to 50-m integrated 
nitrate accurately reproduced observed values. Annual primary production was estimated at -190 gCm' 
yr 1 and varied by no more than ±5% in any year. This estimate is consistent with recent observations but is 
much greater than earlier estimates, indicating that carbon cycling in the North Pacific is much more 
important to the global carbon budget than previously thought. Significant interannual variability in sea 
surface temperature, Ekman upwelling, mixed layer depth, and surface nitrate concentration had little 



impact on productivity. The model results also indicated that the nitrate supply to the euphotic zone is very 
sensitive to Ekman upwelling and that amplification of the wind stress curl results in complete nitrate 
depletion when the winds are persistently downwelling favorable. 

The 30-year simulation at OWS P was repeated using the spectral model (Figure 6) with the zooplankton 
mortality rate, g, adjusted to 0.07/day. This simulation produced improved seasonal chlorophyll a profile 
climatologies (Figure 7) where the elevated values in the boreal summer and fall are more accurately 
reproduced. While there are no profile data available for NH 4 , the seasonal profiles shown in Figure 8 
display subsurface maxima with zero values at depths which are independent of the bottom boundary 
condition. In the simulation of McClain et al. [1996], the subsurface maxima were solely a result of the 
small value assigned as the bottom boundary condition. Finally, Figure 9 illustrates the climatological 
annual time series of 0-60 meter average chlorophyll a (or P), Z, NO,, and NH 4 . The data displayed for 
NO, were taken from Frost [1991] and indicate that the simulated values are slightly high during the winter, 
but are in excellent agreement for the other seasons. 


2.2 Equatorial Pacific Warm Pool at 165° W 

A four-year simulation (1987-1990) of biological processes in the equatorial Pacific Warm Pool was 
performed using the 4-component (phytoplankton, zooplankton, nitrate, and ammonium) spectral ecosystem 
model described above with the input parameter values provided in Table 1. 

The physical parameters ( u , v, w. T) required by the ecosystem model, which are shown in Figures 10-13 
were derived from the OGCM of Murtugudde and Busalacchi [1997], Surface downwelling spectral 
irradiance was estimated using the algorithm described in Section 1.1.3 with ISCCP monthly mean cloud 
fraction data. The vertical diffusivity (Figure 14) was computed from the n, v, and T profiles using the 
method of Pacanowski and Philander [ 1981 J. During 1987—1990, three deep mixing events occurred, all in 
the spring or early summers of 1988-1990. The observed w, v, and T time series are shown in Figures 15. 
16 and 2, respectively. Zonal velocities from the OCGM are in general agreement with observations in 
terms of the magnitudes of the values and the phasing of the flow reversals in the upper 150 m. Meridional 
flow is relatively weak and the model does not emulate the observed values very well. The vertical velocity 
field exhibits a shallow upwelling maxima at about 50 m with alternating periods of upwelling and 
downwelling at depth. Particularly strong upwelling below 100 m occurred in late 1987 and 1990. The 
OCGM temperature field resembles the observed structure but is cooler by 1-3°C with larger differences 
occurring at depth. 

Figure 17 summarizes the ecosystem response to physical forcing. The results clearly indicate that 
upwelling at 100 m is closely associated with surface blooms. For instance, the 100 m upwelling events in 
late 1987 and 1990 produced dramatic increases in the surface layer values of all four ecosystem 
components, whereas the spring-summer deep mixing events (K, time series) do not seem to incur a 
response in any of the ecosystem quantities. 
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Figure 6. Diagram of four-component ecosystem model showing the bio-chemical interactions between 
components. 




Depth (meters) Depth (meters) 


24 


Phyto(mg Chl-a rrr 3 ) Phyto(mg Chl-a rrr 3 ) 

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.1 0.2 0.3 0.4 0.5 0.6 



0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.0 0.1 0.2 0.3 0.4 0.5 0.6 



Model(solid), Observed(dashed), Std Dev.(dotted) 


Figure 7. Seasonally averaged vertical profiles of phytoplankton concentration predicted by the ecosystem 
model at 50°N and 145°W (OWS P). The model profiles are the solid lines, the observed values are 
represented by the dashed lines, and the the standard deviation boundaries are shown by the dotted lines. 
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Figure 8. Seasonally averaged vertical profiles of ammonium concentration predicted by the ecosystem 
model at 50°N and 145°W (OWS P). 
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Figure 9. Composite time series plot of 30-year (1951-1980) vertically averaged (0-60 m) climatological 
concentrations of phytoplankton, zooplankton, nitrate, and ammonium. Note the large seasonal cycle of all 
bio-chemical components. The climatological seasonal cycle of nitrate concentration compares well with 
observed data (blank circles). 












Meridional Velocity in cm/s (OGCM Simulation) 



Figure 11 . Four-year (1987-1990) time series of monthly-averaged OGCM meridional velocity at 165°E 
the equator. 
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Figure 12. Four-year (1987-1990) time series of monthly-averaged OGCM vertical velocity at 165°E on 
the equator. 
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Figure 13. Four-year (1987-1990) time series of monthly-averaged OGCM temperature at 165°E on the 
equator. 
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Figure 14. Four-year (1987-1990) time series of monthly -averaged vertical diffusivity (K v ) at 165°E on the 
equator. The values of K v were calculated using the method of Pacanowski and Philander [3981] and 
physical fields simulated by the OGCM. 
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Figure 15. Four-year (1987-1990) time series of monthly-averaged TOGA-TAO meridional velocity at 
165°E on the equator. 
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Figure 16. Four-year (1987-1990) time series of monthly-averaged TOGA-TAO zonal velocity at 165°E 
on the equator. 
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Figure 17. Composite four-year (1987-1990) time series of vertically averaged (0-50 m) phytoplankton, 
zooplankton, nitrate, and ammonium concentrations, combined with the vertical diffusivity and vertical 
velocity time series at 100 m originating from the OGCM-driven Warm Pool simulation. 
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TABLE 1 

Ecological Model Variables and Input Parameter Definitions and Values. "Original” refers to the values used in 
McClain et al. ( 1996) and "spectral” are values tuned for effects of nitrification and remineralization. 


Symbol OWS P OWS P Warm Pool 

(original) (spectral) 

P 

Z 

NO 3 

NH.i 

G 

7tl 

712 

1 



0.2 

0.2 

0.1 


0.1 

0.1 

0.1 


13.0 

13.0 

0 


0.1 

0.1 

0.1 

m 

0 

0.1 

0.1 

Go 

0.851 

0.851 

0.851 

kgp 

0.0633 

0.0633 

0.0633 

T| 

0.02 

0.02 

0.02 

krp 

0.0633 

0.0633 

0.0633 

Ikma* (<6()m) 

25.0 

25.0 

400.0 

Ikmax (>60m) 

250.0 

160.0 

400.0 

Kn03 

1.0 

1.0 

1.0 

KnH4 

0.1 

0.1 

0.5 

pk 

3.0 

5.0 

5.0 

Smax 

1.0 

1.0 

0 

chla:N 

1.0 

1.0 

1.0 

g 

0.04 

0.1 

0.04 

Rm 

5.0 

4.75 

7.0 

A 

0.8 

0.8 

0.8 

TH 

0.0 

5.0 

0.0 

y 

0.3 

0.3 

0.3 

Tzo 

0.019 

0.019 

(reformulated) 

krz 

0.15 

0.15 

(reformulated) 

az 

0.5 

0.6 

0.6 

a P 

0.5 

0.6 

0.6 

A 11 

A max 


2.0 

2.0 

Dmin 


0.0095 

0.0095 

Kd 


0.036 

0.036 

Cpel 


0.8 

0.8 

KVbot 

17.3 

17.3 

8.65 


Definition 
Derived Quantities 

Phytoplankton nitrogen. mg-at-N/m 
Zooplankton nitrogen. mg-at-N/m 
Nitrate, mg-at-N/m 3 
Ammonium. mg-at-N/m 3 
Phytoplankton growth rate. 1 /day 
Regenerated production fraction 
New production fraction 
Zooplankton grazing rate, 1/day 

Initial Conditions 

Initial P concentration. mg-at-N/m 3 
Initial Z concentration. mg-at-N/m 
Initial surface nitrate concentration, mg-at-N/m* 

Initial ammonium concentration. mg-at-N/m 

Phytoplankton Input Parameters 
Death rate. 1/day 
Growth rate at ()°C, 1/day 
Temperature sensitivity of algal growth. 1/°C 
Respiration coefficient 

Temperature sensitivity of algal respiration. 1/°C 
Maximum photoacclimation parameter, pEin/nr/s 

Half saturation for nitrate uptake, mg-at-N/m 3 
Half saturation for ammonium up-take, mg-at-N/m* 

Ammonium inhibition of nitrate uptake, dimensionless 
Maximum sinking rate, m/day 
Chlorophyll-a:Nitrogen ratio, mg/mg-ut-N 

Zooplankton Input Parameters 
Death rate, 1/day 
Maximum grazing rate, 1/day 
Ivlev constant, m 3 /mg-at N 
Minimum C threshold for grazing, mgC/m 
Unassimilated ingested ration, dimensionless 
Respiration rate at 0°C, 1/day 
Temperature sensitivity of respiration. 1/°C 

Nutrient Input Parameters 
Fraction of dead zooplankton converted to ammonium 
Fraction of dead phytoplankton converted to ammonium 
Maximum rate of nitrification, nmo 1/day 
Minimum irradiance (300-470 nm) for nitrification, W/m" 
Half-saturation dosage for nitrification photoinhibition. W/nT/nm 
Fecal pellet remineralization fraction 

Physical Parameters 

Minimum (bottom) eddy diffusion coefficient, nr /day 
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APPENDIX A 

The Interactive Data Analysis Package (IDAPAK) 

IDAPAK is a software package designed at NASA GSFC to perform user-specific functions in connection 
with the preprocessing and analyses of geophysical data sets. IDAPAK was originally developed tor 
personal computers (PCs) to complement PC-SEAPAK (Firestone et al., [1989]; Firestone et al., 1990; 
Darzi et al.. [ 1991]), which did not support meteorological and hydrographic data analysis. Over the past 
few years, PC-IDAPAK has been ported to UNIX (Tai and McClain [ 19931) where its functionality has 
been expanded by implementing the package under the commercially available Interactive Data Language 
(IDL) environment by Research Systems, Inc. 


A.l Hardware and Software Configurations 

While most of IDAPAK programs are coded in FORTRAN-77 and C, the package makes extensive use of 
IDL, along with its graphical user interface (GUI), which is available for most UNIX systems. IDL is a 
powerful interpreter computer language which greatly expands the data analysis and display capabilities of 
IDAPAK. It provides many graphical and image display functions and also incorporates many data 
transformation and statistical application programs for image processing, data visualization, and scientific 
data analysis. The package can access data sets from a variety of sources and convert them into a simple 
columnwise ASCII format which is compatible with the other functions in the package and other PC and 
UNIX applications. All the input and output files in connection with ECOID are in IDAPAK format. 


A.2 Functions 


The modules in IDAPAK are grouped into six different categories according to the characteristics of each 
module or program. They are: 


• ACCESS (Data Access) 

• PROCESS (Multivariable File Manipulation) 

• ANALYSIS (Statistical Analysis and Derived Product Generation) 

• MODEL (Oceanic Ecosystem and Radiative Transfer Models) 

• DISPLAY (Data Display) 

• ANIMATION (Data Animation) 


Table Al shows the current implementation of IDAPAK. Each module is an independent function. Most 
sources of meteorological and oceanographic data use different formats making the usage of the data and 
the development of application software more difficult. Olsen and McClain [ 1992] report on one effort to 
gather key oceanic data sets and archive them in a common format (the NASA Common Data Fonnat, 
CDF) at the NASA Climate Data System (NCDS). This activity paralleled the development of the 
SEAPAK (McClain et al. [ 1992]) environmental applications module and interdisciplinary research using a 
variety of data sets (McClain et al., [1990]; Brock et al„ [1991]; McClain and Firestone. [1993]) most of 
which were in CDF. NCDS has been phased out and support for many of the ocean data sets discontinued 
as NASA is now focused on the development of the EOS Distributed Active Archive Centers (DAACs). 
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Thus, most data sets presently being obtained are no longer in CDF, but are in the original formats. The 
programs under ACCESS reflect some of the more recent data sets used by the authors for ocean research. 

The following subsections provide descriptions of each IDAPAK module. 

A.2.1 ACCESS 

The access programs decode data and convert them into an ASCII format designed for IDAPAK This 
ASCII format arranges data in column-row form with a small information header describing the number of 
columns, the format of each column and the content (variable descriptor) of each column. Data files with 
this format are transportable and easily edited. For example, they can be downloaded to a PC and 
processed using any other analysis or graphics software. The present suite of conversion programs can 
handle hydrographic station data from NOAA’s NODC (National Oceanographic Data Center), data in 
netCDF (network Common Data Format), data from TOGA-TAO array, data generated by FOXPRO (a 
Microsoft database product) and several other data sets associated with specific field experiments or 
investigators. TOGA is an ongoing project for investigating the interactions between the tropical 
atmosphere and ocean. The TAO array of buoys are deployed across the tropical Pacific and some include 
solar irradiance measurements. TOGA-TAO data are either in ASCII or netCDF format. Once read by the 
ACCESS programs, the data can be processed, analyzed, and displayed using the other modules in 
IDAPAK. Table A2 describes each specific ACCESS function. 

A.2.2 PROCESS 

The PROCESS programs perform basic data and file manipulation operations such as copying, deleting, 
renaming, editing, cutting and pasting, listing, printing, merging, or appending different files into one file, 
reformatting, transposing, and chronological sorting. Table A3 describes each specific PROCESS function. 

A.2.3 ANALYSIS 

The programs tor basic statistical analysis such as auto- and cross-correlation, linear regression, fast Fourier 
transformation, and maximum entropy spectral analysis are collected under the ANALYSIS module. There 
are also functions for smoothing, interpolation using different algorithms (linear and cubic spline), and 
simple columnwise algebraic operations for multi-variable files. Table A4 describes each specific 
ANALYSIS function. 

A .2.4 MODEL 

The MODEL module incorporates a variety of analytical and numerical model products. DERIVED V is 
a program that computes 20 meteorological and oceanographic derived quantities from basic input (see 
Table A5). 

SOLAR G&C, a solar irradiance model by Gregg and Carder [1990], estimates the clear sky surface 
spectral downwelling solar irradiance at 1 nm resolution from 350 to 700 nm. SOLAR F (Frouin et al., 
[1989]), estimates broad-band clear sky solar irradiance for 250-4000 nm, 350-700 nm and 400-700 nm. 
SOLAR_G computes the total clear-sky solar irradiance (Garwood, personal communication). 

COUPLED_M includes 3 one-dimensional marine ecosystem models which differ mostly in the way that 
the physical parameters are calculated. The first one couples mixed layer dynamics [Garwood, 1977], 
biological processes [McClain et al., 1996] and the surface irradiance models (Frouin et al., [1989]; Gregg 
and Carder, [1990]) . The model is driven by surface meteorological variables (wind speed, cloud cover, 
air temperature, and dew point temperature). These variables are used to determine the surface latent, 
sensible, short wave and long wave heat fluxes, and the wind stress. The subsurface light field is 
determined using the vertical profile of chlorophyll as estimated from the biological model. 
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Table Al. IDAPAK programs listed by program category. 


ACCESS 

PROCESS 

ANALYSIS 

ARRIGO 

ADDJUL1AN 

AREAVGIDA 

BISHOP_PAR 

: APPEND 

AVERAGING 

BRS_CDF 

BOUND 

AVG.FILES 

BRS_DAAC 

CONV_FMT 

BANDPASS 

BRS_HDF 

COPY 

BINNING 

BRS.NETCDF 

DEL_ROW 

COLUMNOPER 

CDFLST2IDA 

EDIT 

CORRELATION 

CREATE 

IMAGE 

FFTIDA 

ENV_QUERY 

INSERT.COL 

FUNCSIDA 

EXCEL2IDA 

LIST 

GETIDAMLD 

FOXPRO 

MERGE 

INTEGRATION 

FRENCH 

RANGE 

INTERPOLA 

GARWOOD 

REFORMAT 

MEMSPCIDA 

GETCDF 

REMOVE 

REGRESSION 

GETDAAC 

RENAME 

SMOOTH 

GETIDLBIN 

RM_9999 

SPLINE 

GETISCCP 

PAGECUT 

STDEV 

GETNODC 

PRINT 

SUMCOL 

GETRAGU 

GETSALT 

GETSUE 

GETWOA 

TAO_RAD 

TAO_RAD_20 

TOGA_TAO 

TIMENV21DA 

SAMPLE 

SELECT 

SORTING 

STRIPHDR 

TRANSPOSE 

SUMIDA 


MODEL 

DISPLAY 

ANIMATION 

DERIVED_V 

SET_WINDOW 

NODC_PROF 

COUPLED_M 

INSIGHT 

IDA_LINE 

SOLARJ3&C 

CREATE.MAP 

IDA_PROF 

SOLAR_F 

VIEW_CZCS 

CZCS 

SOLAR.G 

VIEW_HDF8 

HDF8 


VIEW_HDF24 

HDF24 


VIEW_SFC 

DAAC.BIN ary 


VIEW_PROF 

VIEW_LINES 

PLOT_CONTRS 

PLOT.SFC 

PLOT_LINE 

PLOT_PROF 

PLOT.SCATTR 

SURF_CONT 

TRACKING 

TOGA_CROSS 

VECT OR_PLOT 

SUE_BINARY 
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Table A2. Description of IDAPAK ACCESS functions. 


ARRIGO 



BRS HDF 


CDFLST2IDA 


CREATE 


env_query 


EXCEL2IDA 


FOXPRO 


GARWOOD 


FRENCH 


GETCDF 


GETDAAC 


GETIDLBIN 


GETISCCP 


GETNODC 


GETRAGU 


GETSALT 


GETSUE 


GETWOA 


TAO RAD 20 


TIMEN2IDA 


Data Formatted by K. Arrigo (NASA/GSFC) 


Data Archived by J. Bishop in HDF SD format 


NASA CDF browse data 


DAAC binary browse data 


netCDF browse data 


HDF browse data or images, and entire array dum 


Converts CDFLIST ASCII data to IDA format 


Creates data using NEDIT 


Inquires CDF data sets stored in OCF at NASA/GSFC 


Converts from Microsoft Excel CSV ASCII data to IDA format 


Extracts data generated by FOXPRO 


Data formatted by R. Garwood 


Reads data generated by French research ships 


Extracts data from two-dimensional CDF ASCII data sets 


Extracts data from DAAC binary files or IDL binary files 


Extracts IDL binary data and converts to 2D CDF ASCII or row/column format 


Converts NASA/Langley ISCCP HDF data to CDF 2D ASCII or row-column formats 


Extracts NODC data from optical disk or CD-ROM 


Extracts OGCM data generated by R. Murtuggude 


Extracts salinity profile data and converts it to ASCII column format 


Extracts SUE surface flux data 


Extracts CD-ROM data archived by WOA 


Extracts TOGA-TAO meteorological data formatted in netCDF or ASCII 


Extracts irradiance data archived from TAO buovs 


Averages 20-minute irradiance data from TAO buoys to hourly data 


Converts SEAPAK TIMENV ASCII data to IDL 
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Table A3. Description of IDAPAK PROCESS functions. 


ADDJULIAN 

Adds Julian date to an IDAPAK file 

APPEND 

Appends IDAPAK files into a single file 

BOUND 

Allows user to specify high and low acceptance bounds for variables in the data file 

CONV FMT 

Converts IDAPAK, CZCS, and CDF files into other formats (ASCII, CZCS. and CDF) 

COPY 

Copies one IDAPAK file into another 

DEL ROW 

Deletes rows from IDAPAK files (start/stop dates must be specified) 

EDIT 

Allows user to edit files on an IDAPAK window 

IMAGE 

Image processor (composites, convertions,histograms, median and filters) 

INSERT COL 

Inserts columns into an IDAPAK file 

LIST 

Lists files on an IDAPAK window (general, IDAPAK, and NODC formats) 

MERGE 

Merges the data columns of IDAPAK and ASCII files into a single tile 

RANGE 

Disabled 

REFORMAT 

Converts IDAPAK files from floating to exponential formats, and vice-versa 

REMOVE 

Deletes a file 

RENAME 

Renames a file 

RM_9999 

Remove record nulls from a file 

PAGECUT 

Allows user to extract a subset of variables from a file 

PRINT 

Prints a file 

SAMPLE 

Sub-samples an IDAPAK or NODC file according to time and geographic bounds 

SELECT 

Selects data from a file based on station ID, date and hour, and prints to another file 

SORTING 

Sorts IDAPAK time series and profile data into yearly, monthly, and seasonal time bins 

STRIPHDR 

Strips the header information from IDAPAK or ASCII files 

TRANSPOSE 

Transposes rows and columns in a file and prints results to another tile 


Table A4. Description of IDAPAK ANALYSIS functions. 


AREAVGIDA 

Computes area-averaged values from a file given geographical boundaries 

AVERAGING 

Computes total,yearly,monthly,daily,12-,6-,3-hourly, and hourly time averages 

AVG FILES 

Combines any number of specified files and averages them into a single tile 

BANDPASS 

Applies a bandpass filter to variables in a file given the lower and high frequencies 

BINNING 

Bins the data in a file according to monthly, seasonal or total means 

COLUMNOPER 

Applies user-specified operators +, -, x, and / to individual columns in a tile 

CORRELATION 

Calculates the cross-correlation between variables within a fde or in separate tiles 

FFTIDA 

Applies a Fast Fourier Transform (FFT) to specific variables in a file 

FUNCSIDA 

Applies arithmetic, logarithmic, and trigonometric operators to columns in a tile 

GETIDAMLD 

Gets mixed-layer depth from a temp record using Levitus or buoyancy methods 

INTEGRATION 

Reads vertical profiles and averages, integrates, and picks values from given levels 

INTERPOLA 

Interpolates time series and profiles at given intervals of time or depth 

MEMSPCIDA 

Performs spectral analysis given filter and bandwidth specifications 

REGRESSION 

Linear regression between two variables within a file or separate files 

SMOOTH 

Filters specific data sets within a file given filter characteristics 

SPLINE 

Fits a spline into vertical profiles or time series given the spline and data specs 

STDEV 

Calculates max, min, sum, mean, and std deviation of all columns in a file 

SUMCOL 

Sums and multiplies data columns in a file given a specific operator 

SUMIDA 

Time series sum (as in AVERAGING, except that it is not divided by no. of points) 
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Table A5. List of quantities that can be computed in IDAPAK given the required input fields. 


1. Water density 

22. Longwave radiation (Garwood for clear skv) 

2. Dew point temperature 

23. Longwave radiation (Berliand and Berliand, 
[1952]) 

3. Vapor pressure 

24. Longwave radiation (Budyko, f 1963]) 

4. Relative humidity 

25. Longwave radiation (Budyko. [ 19631 ((Clear skv) 

5. Mixing ratio 

26. Surface irradiance (Budyko, [ 19631) 

6. Potential temperature 

27. Surface irradiance (Reed, f 19771) 

7. Equivalent potential temperature 

28. Surface irradiance (Dobson & Smith, [1988] ) 

8. Zonal & meridional wind components 

29. Clear sky surface irradiance (ocean) (Davis) 

9. Zonal & meridional wind stress ( x , y ) 

30. Clear sky surface irradiance (ocean) (Atwater) 

10. Total wind stress (t) 

31. Clear sky surface irradiance (land)(Davis) 

ii. 

32. Clear sky surface irradiance (land) (Atwater) 

12. Zonal & meridional Ekman transport 

33. Visibility (km) 

13. Total Ekman transport 

34. Water content (cm) 

14. Ekman depth 

35. Albedo 

15. Sensible heat flux (SH) 

36. Solar zenith angle 

16. Latent heat flux (LH) (4 inputs) 

37. Cosine of zenith angle 

17. Latent heat flux (LH) (5 inputs) 

38. Solar altitude 

18. Bulk coefficient for SH 

39. Julian date 

19. Bulk coefficient for LE 

40. Clear sky surface irradiance (Ocean) (Frouin) 

20. Longwave radiation (Garwood, 3 input 
values) 

41. Clear sky surface irradiance (Land) (Frouin) 

21. Longwave radiation (Garwood, 4 input 
values) 



The second model is similar to the first model, except that the mixed layer dynamics are simulated by the 
Chen et al. [1994] one -dimensional model. As mentioned in Sections 1.3 and 1.4, the mixed layer model of 
Chen et al. [1994] was modified by adding flux and diffusion adjustment terms for preventing excess 
heating at the bottom of the model domain, which propagates upward to the surface. Besides this 
modification, biological processes are reformulated to include nitrification of ammonium and 
remineralization of fecal pellet production. The vertical velocity field is either derived from the ocean 
model (Murtugudde et al. [1996]) or calculated from an empirical formula (McClain and Firestone [1993]). 
The algorithm for the estimation of diffusivity in the water column uses the method of Pacanowski and 
Philander [1981], This method requires current velocity and potential temperature. The TAO buoys 
deployment over the tropical Pacific provide daily records of these variables; however, the temporal data 
set coverage is very limited. To have a longer period of simulation, the model was upgraded into a third 
version that allows the ingestion of profiles of current, temperature, and vertical velocity computed by the 
OGCM (Murtuggude et al., [1994]). All these models generate time series and profiles of different 
variables for diagnostics and analysis purposes. Tables A6, A7, A8, A9, and A10 show examples of model 
output stored on five different output files that can be viewed and analyzed with the IDAPAK Graphic User 
Interface (GUI) functions described in Section A.2. These files can be easily modified in the corresponding 
FORTRAN modules to satisfy specific user needs and future model upgrades. 
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Table A6. Time series of diagnostic and predicted variables (file fc namef.asc’ in Section B2.2) 


1. Mixed layer depth 

14. Zooplankton 

2. Predicted sea surface temperature 

15. Ammonium 

3. Observed sea surface temperature 

16. Nitrate 

4. Air temperature 

17. Radiation at surface 

5. Salinity 

18. Radiation at mixed layer 

6. Q (sensible heat) 

19. Averaged mixed layer depth 

7. Q (longwave radiation) 

20. Gross production 

8. Q (latent heat) 

21. Net production 

9. Q (radiation just below surface) 

22. New production 

10.Q (radiation of clear sky) 

23. F ratio 

1 1 . Q (radiation at surface) 

24. Radiation at mixed layer of OGCM 

12. Radiation at mixed layer depth 

25. Mixed layer depth of OGCM 

13. Phytoplankton 

26. PAR at surface 


Table A7. List of variables in the output file containing vertical profiles generated by the coupled 
model (file ‘namep.asc’ in Section B2.2). 


1. Depth 

13. Vertical Diffusivity 

2. Julian date 

14. Phytoplankton Respiration rate 

3. Temperature 

15. Vertical velocity 

4. Salinity 

16. C:Chla 

5. Phytoplankton 

17. Zooplankton respiration rate 

6. Zooplankton 

18. 24-hour running mean of subsurface light (E ) 

7. Ammonium 

19. Llim 

8. Nitrate 

20. Nlim 

9. 24-hour averaged growth rate 

21. Beta term 

10. Grazing rate 

22. NOjJim 

ILF ratio 

23. NH 4 _lim 

12. Profile of PAR 

24. Integrated daily gross production 
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Table A8. List of variables in the output file containing the time series of integrated data (file 
‘namet.asc’ in Section B2.2). 


1 . Phytoplanton integrated to 50 m 

1 6. F ratio for column 

2. Zooplankton integrated to 50 m 

17. Sensible heat flux 

3. Ammonium integrated to 50 m 

18. Heat flux due to back radiation (longwave) 

4. Nitrate integrated to 50 m 

19. Latent heat flux 

5. Phytoplankton integrated for column 

20. Clear sky total radiation just above the surf. 

6. Zooplankton integrated for column 

21. Cloudy sky total radiation just above the surf. 

7. Ammonium integrated for column 

22. Cloudy sky total radiation just below the surf 

8. Nitrate integrated for column 

23. Gregg-Carder PAR just below the surface 

9. Gross production for 50 m 

24. Gregg-Carder PAR at the mixed-layer depth 

10. Net production for 50 m 

25. Frouin clear-sky PAR just above the surface 

1 1. New production for 50 m 

26. Frouin clear-sky PAR just below the surface 

12. Gross production for column 

27. Frouin cloudy-sky PAR just above the surface 

13. Net production for column 

28. Frouin cloudy-sky PAR just below the surface 

14. New production for column 

29. Dobson clear sky radiation 

15. F ratio for 50 m 



Table A9. More time series of predicted and diagnostic model variables (file ‘namen.asc’ in Section 
B2.2). 


1 . Mass loss due to grazing 

10. Observed temperature at 250 m 

2. Column-integrated pet of Phyto transformed into ammonium 

11. Predicted temperature at 250 m 

3. Column-integrated pet of Zoo transformed into ammonium 

12. Ratio of light penetration to total 
surface irradiance (y) 

4. Nitrate concentration at the bottom 

13. Avg growth within euphotic zone 

5. Vertical advection of N0 3 at the bottom 

14. Avg Phyto within euphotic zone 

6. Sinking rate of Phyto at the bottom 

15. Avg NO} within the top 5 depths 

7. Total mass at each time step (Phyto+Zoo-t- NOi+NHj) 

16 Avg N0 3 between 56-65 m. 

8. Total mass change at each time step 

17. Frouin lo Gregg-Carder PAR ratio 

9. Total mass loss at the bottom 



All the models and formulas need initial input data of various kinds. This information is passed to the 
model or formulas through the GUI. The input data and results from these models are all preprocessed 
using other IDAPAK modules. Post-processing and display of model outputs are also accomplished using 
IDAPAK modules. The user may access data from different data sources and generate the inputs required 
by the models using programs in PROCESS and ANALYSIS modules. The" results of models may be 
displayed or animated by functions in DISPLAY and ANIMATION. 

A.2.5 DISPLAY and ANIMATION 

DISPLAY and ANIMATION have functions for displaying and animating the column— row formatted data 
and image files. Since IDL provides a very comprehensive graphical display and image processing library, 
all the programs in these two modules are coded in IDL. The display module allows the user to display 
one-dimensional time series data, contour two-dimensional fields, and annotate graphs interactively. For 
irregularly distributed horizontal data, the user can generate contours on a variety of map projections. All 
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DISPLAY programs allow the user to generate PostScript™ output. The ANIMATION programs support 
looping ASCII profile data and CZCS image files. 


Table A10. Diagnostic file containing daily profiles of all phytoplankton, zooplankton, nitrate and 
ammonium terms in (41) through (44) (file ‘diag bal name.dat 1 in Section B2.2). 


1. Depth 

2. Phyto time change (3P/3t) 

3. Phyto vertical advection (w3P/dz) 

4. Phyto sinking rate (3SP/3z) 

5. Phyto vertical diffusion (k v d 2 P/dz") 

6 . Phyto growth rate (GP) 

7. Phyto mortality rate (mP) 

8 . Phyto respiration rate (r p P) 

9. Zoo grazing on Phyto (IZ) 

10. Zoo time change (dZ/dt) 

1 L Zoo vertical advection (wdZ/3z) 

12. Zoo vertical diffusion (k v 0 2 Z/0z~) 

13. Zoo growth ([ 1-ylIZ) 

14. Zoo death rate (gZ) 

15. Zoo respiration rate (r z Z) 


16. NQ 3 time change (ONQj/dt) 

17. NO 3 vertical advection (w3 NQVOz) 

18. NO 3 vertical diffusion (k v D 2 NOydz 2 ) 

19. NOj uptake by phytoplankton (k^GP) 

20. Ammonium nitrification (A n ) 

2 L NH 4 time change (3NH 4 /0t) 

22. NH 4 vertical advection (wd NH 4 /dz) 

23. NH 4 vertical diffusion (k v 3~ NH 4 /3z 2 ) 

24. NH 4 production from phyto mortality (a D mP) 

25 NH 4 p roduction from phyto respiration (r p P) 

26. NH 4 loss due to phyto growth (~n { GP) 

27. NH 4 production from zoo mortality (a z gZ) 

28 . NH 4 production from zoo respiration (r 7 Z) 

29. NH 4 production from zoo fecal pellets (CpeiylZ) 

30 NH 4 loss due to nitrification (-A n ) 
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APPENDIX B 

Model Configuration and FORTRAN Routine Modules Used by the 
One-Dimensional Ecological Model (ECOID) 

This appendix provides descriptions of the model configuration and FORTRAN routine modules used by 
two versions of the ecosystem model. Version 1 consists of a coupled Garwood [1977] MLM with the 
spectral version of the 4-component ecosystem model. Version 2 uses the same 4-component ecosystem 
model but the physical forcing fields (u,v,\\\T) are either derived from observations (e.g., TOGA-TAO 
array) or from an independently run OGCM. Since both versions share the same bio-optical algorithms, a 
certain degree of overlap in the descriptions of the two model versions is inevitable. The following sections 
address details of model configuration and usage. 

B,1 ECOID with Garwood MLM Coupling 

This section describes the input parameters, nominal values, input and output files, program flowchart, and 
FORTRAN routine names and functions for the ecosystem model using the Garwood [ 1977 1 MLM. 

B.1.1 Input Parameters and Nominal Values 

The input parameters and nominal values used by the ecosystem model are read in from a file (filename. par) 
specified in the command line when executing the model run: 

c k ecold < filename. par & 

The initial and forcing fields required for the multi-year runs are read in from other files specified in 
Section B1.2. The format for all parameters in ‘filename. par’ is provided in the main FORTRAN module 
‘opblpc.f . The following table summarizes each of the input parameters contained in the file. 

IC= Types of initial temperature (T) and salinity (5) profiles 

-1: use default T{ z) and S(z) 

0: use analytic T{ z) and default S(z) 

1 : use Garwood T(z) and default S(z) 

2: user’s provided 7\z) and default S{ z) 

3: user’s provided T{ z) and S( z) 

4: use T( z) and S( z) produced in last run 

IB= Boundary condition source 

0: analytical forcing test case for MLM 
1: read boundary values from external sources 

IE= Upwelling option: 

0: use constant vertical velocity We 
1 : evoke Ekman upwelling calculation 


CB= Biological model option: 

0: no biological model evoked 

1: couple mixed layer model with biological model 

FB= 


Switch to enable bio-absorption of light 
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IP= 


DD= 


DB= 


CV= 


HR = 24 
HP = 24 
v= 2.79 
uo= 0.34 


a_type= 


Cldbb=0.5() 
Cldaa=0.53 
shift=5.0 
dt=0.3600D+04 
dz =0.1000D+0I 
phs= 0.0 
per= 365 
A75= 0 


NST= 


AT= 0.02 
AB= 0.02 


0: no 
1: yes 

Switch for initial nutrient concentration file: 

1: use default values 

2: read initial nitrate and ammonium concentrations from units 7,8 
3: read initial nitrate and ammonium concentrations from hot start file (7) 

4: read initial nitrate concentration from unit 7 and use default ammonium 
values 

Switch to print results: 

0: do not print results 
1 : print results 

Switch to print daily-averaged time series of surface and depth-averaged values: 
0: do not print time series 
1 : print time series 

Switch to control writing of hot start file: 

0: do not write hot start file 
1: write hot start file 

Controls frequency (in hours) of output of model variables (except for profiles) 
Controls frequency (in hours) of output of vertical profiles 
Atmospheric water content (cm) 

Ozone concentration (Doobson) 

Aerosol type: 

0: Continental 
1: Oceanic 

Cloud correction coefficient 
Cloud correction coefficient 

Shift angle (degrees) for the periodic seasonal cloud effect 
Model time step (seconds) 

Vertical grid resolution (meters) 

Phase for analytical vertical motion parameters 
Period for analytical vertical motion parameters 
Internal wave amplitude at 75 meters 

Controls equation of state 

0: constant thermal expansion coefficient 

1: uses temperature-dependent thermal expansion coefficient 

Diffusion coefficient at the top of the mixed layer (in cm 2 fs) 

Diffusion coefficient at the bottom of the mixed layer (in cm 2 /s) 


NV= Switch to control diffusion in MLM 

1 : temperature diffusion only 
4: diffusion of all variables 

H=350.0 Maximum depth for model calculations 

Pl= 1 Coefficient for ML entrainment equation 



P2= 1 
P3= 1 
M3 = 0.5 
K= 0.004 
AKD=0.0525 
ZONE=25.0 
BOTKV=8.65 
Mlcr= 1 
Mlv= 0.1 
Mcr= 30 
We= 0.25 
LAT= 50.0 
LON= -145.0 

IW= 


m= 0.1 
G0= .8511 
Kno3= 1.0000 
Knh4= .5000 
c= 0.02 
Ideep = 25.0 
Ishal= 160.0 
mixed= 60.0 
NH4= 0. 1 
N03= 13.0 
P= 0.2 
Z= 0.1 
pk= 5.0 
Knir= 2.5 
K v = 1.0 
CHLN= 1.0 
DI = 2 

cs= 2 
fc=1.0 
Smax= 1 .0 

g= 0.1 

Rm= 4.75 
LAMBDA= 0.8 
GAMMA= 0.3 
rz0= 0.0190 
krz= 0. 1 5 
Znh4 = 0.6 
Pnh4 = 0.6 
TH = 5.0 
MOVING= 24 
START= 19510101 
NDAYS= 10957 

IMLD- 1 


Coefficient for ML entrainment equation 
Coefficient for ML entrainment equation 
Coefficient for heat flux equation 
Extinction coefficient 

Coefficient for the vertical profile of diffusion used by the bio model 
Depth of transition layer (meters) 

Minimum vertical diffusivity (Kv in nr/day) 

Controls number of variables to be printed (set to 12) 

Criterion coefficient to find mixed-layer depth 
Minimum} mixed-layer depth in meters 
Upwelling velocity (m/day) 

Latitude in degrees (North is positive) 

Longitude in degrees (East is positive) 

Control for upwelling velocity: 

1: constant We 

2: read in We from external source 

Phytoplankton death rate, d 1 
Phytoplankton growth rate at 0°C, d 1 
Half saturation of nitrate uptake, mg-at N m' 3 
Half saturation of ammonium uptake, mg-at N m 
Ratio of phytoplankton respiration to growth 

Maximum photoacclimation parameter (> 60 m), pmol quanta/nr/s 
Maximum photoacclimation parameter (< 60 m), pmol quanta/nr/s 
Critical depth (meters) for Ishal and Ideep calculation 
Default initial NH4 value (if input file not assigned) 

Default initial N03 value (if input file not assigned) 

Initial Phytoplankton concentration, mg-at N m 
Initial Zooplankton concentration, mg-at N m 
Nitrate-sensitive coefficient due to ammonium 
Attenuation coefficient for Near Infrared (m 1 ) 

Minimum value for K v (m'/d) 

Chlorophyll-a:Nitrogen ratio, mg/mg-at N 

Index to control the method to compute surface diffusivity 

Index for choosing the method to compute the diffusivity profile 

Multiplicative factor to apply to vertical diffusivity 

Sinking speed (m/d) 

Zooplankton death rate, d 1 
Zooplankton maximum grazing rate, d 1 
Ivlev constant, m 3 /mg-at N 

Unassimilated zooplankton ingested ratio, dimensionless 
Not used - reformulated 

Temperature sensitivity of zooplankton respiration, °C 1 
Fraction of dead zooplankton converted to ammonium 
Fraction of dead phytoplankton converted to ammonium 
Minimum C threshold for zooplankton grazing, mgC/m 3 
Moving average window in hours for light energy 
Starting date for calculations, year/month/day 
Total number of days of simulation 

Switch for depth integration limits (T1 ) 

1: surface to mixed-layer depth 
2: surface to bottom 
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T2= 50.0 
T3= 350.0 
DLEVEL= 350.0 
INTV= 10.0 
BLKCH=0.0015 
BLKCE=0.0015 
NDG=2 
NFX=1 


Integration depth limits (meters) 

Integration depth limits (meters) 

Maximum depth level for model output (meters) 

Depth interval for output (meters) 

Transfer coefficient for sensible heat flux 
Transfer coefficient for latent heat flux 
Controls calculation of drag coefficient for wind stress 
Controls calculation of heat flux transfer coefficient 


NLW=1 


NH4OX=2.0 

HALF=0.036 

ZtoNH4=0.8 

CLDF=1.0 

CLC=0.001289 


Switch to control longwave radiation scheme 
1: Tseng method (W/nr) 

2: Budyko method (Kcal/cm 2 /month) 

Maximum rate of nitrification (nmol/day) 

Half-saturation dosage for nitrification photoinhibition (W/m 2 /nm) 
Remineralization coefficient for fecal pellets 
Multiplicative factor to apply to cloud cover 
Cloud cover coefficient 


DRAG= 1 


YR=195 1 
JULIAN=1 


Controls drag coefficient scheme 
1 : Amorocho and Devries 
2: Constant coefficients 

Beginning year 
Beginning Julian day 


B.1.2 Input and Output Files 

The following logic units (LU) and input/output file names are specified in ‘filename. par’ after the last 
parameter input described in section Bl.l. The file names and data sets that they contain are only given as 
examples; other data sets for different time periods and geographical areas can be specified by the user. 
Note that LU=TMP means temporary logic unit. 

Logic units (LU) for input/ output files 

1 2 3 11 12 13 14 7 8 

Input files 

/user_pathname/xbt 250.dat LU=0 1 ; Initial temperature profiles from XBTs 

/user_pathname/papa5 181.dat LU=02; Environmental data from OWSP ( 195 1-1981) 

/user_pathname/ekmanl.dat LU=03; Pre -calculated Ekman vertical velocities (w e ) 


Output Files 

/user_pathname/namef.asc LU=1 1; Daily time series of selected diagnostic and predicted 

model variables (see Table A6) 
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/user_pathname/namep.asc 

LU=12; 

Daily time series of vertical profiles for selected model 
variables (see Table A7) 

/user_pathname/namet.asc 

LU=13; 

Daily time series of column-averaged and surface 
values for selected model variables (see Table A 8 ) 

/user_pathname/namen.asc 

II 

More daily time series of selected diagnostic and 
predicted model variables (see Table A9) 

Input files 



/user_pathname/decno350.dat 

LU=07; 

Climatological nitrate (NO 3 ) concentration profile 

/user_pathname/winnh4.dat 

LU=08; 

Climatological ammonium (NH 4 ) concentration profile 


Output Files 

/user_pathname/lastdyn.bin LU=TMP; Binary file containing dynamic variables from last time 

step to be used for hot starting the model 

/user_pathname/lastbio.bin LU=TMP; Binary file containing biological variables from last time 

step to be used for hot starting the model 


Input files 

/user_pathname/uvatmodd.dat LU=TMP; Spectral absorption coefficients for Gregg-Carder model 

/user_pathname/wv_chl_absp.dat LU=TMP; Attenuation coefficients due to chlorophyll 
/user_pathname/asnh4.dat LU=TMP; Ammonium nitrification coefficients 

B.1.3 Ecosystem Model Flowchart 

Figure Bl.l shows a flowchart for the first version of the ecosystem model, which consists of a one- 
dimensional MLM (Garwood, [1977]) coupled with the 4-component biological model. The flowchart 
identifies each main FORTRAN module in a separate box, with the arrows indicating input/output data 
flows and time stepping loops. 
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CALL RETR 

> COMPUTES ML 
RETREAT 


* 

CALL SHALLOW 

LIMIT MINIMUM 
MLD TO 101 CM 


CALL RMODE 

VERIFY THAT ML 
RETREAT 
CONSERVES B,T 


▼ 

CALL ABSRP1 

ESTIMATE 
TOTAL ENERGY 
AFTER 

ENTRAINMENT/ 

DETRAINMENT 


* 

CALL 

GREGG_CARDER 

COMPUTE PAR AT 
EACH DEPTH 


CALL 

WATERRAD 

COMPUTE 
WATER COLUMN 
RADIANCE 




CALL 

BIOABSORB 

COMPUTE BIO 
FEEDBACK 


t 

CALL EKMW 

COMPUTE 

EKMAN 

UPWELLING 
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Figure Bl.l. Flowchart of ecosystem model coupled with the Garwood MLM. 
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B.1.4 Routine Names and Functions 

To better describe the functions of each FORTRAN module in the flowchart of Figure Bl.l, a short 
explanation of each FORTRAN module is provided in Table Bl.l. The list of modules in Table Bl.l is 
more extensive than the modules appearing in Figure Bl.l because some FORTRAN modules contain calls 
to subordinate FORTRAN routines and functions. The executable file for the ecosystem model is built by 
compiling and linking all relevant FORTRAN modules using a ‘Makefile’ provided with the source code 
package. 

Table Bl.l. Names and functionality of ecosystem model FORTRAN routines. 


Routine Name 


Description 


ABSRP1 


Absorption of radiation below the surface 


ALGO 


Contains: 

RIGHTHANDPZ,RIGHTHANDNN,PLEFTNEWPZ,PLEFTNEWNN, GAUSS 


BSMN 


Computes critical values of buoyancy to be used in entrainment calculation 


BIOABSORB 


Computes light absorption due to biology 


ATMODD 


Model of atmospheric transmittance of solar irradiance through marine atmosphere 


COMP ESTAR 


Computes running daily mean of light at each depth (E ) 


COMPUTE KV 


Computes vertical eddy diffusivity coefficient for biological model equations 


DAYTIM 


Converts Gregorian to Julian, and Julian to Gregorian dates 


DIFF 


Computes vertical diffusion term (see also GARDIFF) 


DIFFUSE 


Name of FORTRAN file containing routine DIFF 


EKMW 


Estimates the Ekman upwelling profile given wind speed and latitude 


ENERGY 


Computes energy balance within the mixed layer 


ENTR 


Computes mixed-layer growth and adjusts profiles of temperature, salinity and buoyancy 


FROUIN 


Computes clear sky solar irradiance using Frouin et al. method 


GAUSS 


Gauss elimination implicit matrix solver for ecosystem model (Crank-Nicholson scheme) 


GARDIFF 


This routine is also part of the vertical diffusion computation 


GREGG 


Computes spectral solar irradiance using the Gregg and Carder method. Contains: 
GREGG CARDER 


GROWTH 


Computes respiration and growth. Contains: 
GROWTH 1 ,GROWTH2,RESIZ 


HFLUX 


Computes the surface heat flux 


ICOND 


Reads initial profiles of temperature and salinity 


IMPMIX 


Adjusts vertical profiles according to the depth change of the mixed layer 


INITPZN 


Initializes Phyto, Zoo, NH 4 , and NQ 3 concentrations 


INIT 


Initializes temperature and salinity fields 


JPMIX 


Relieves gradient Richardson number instability using J. Price’s criterion 


OPBLPC 


Main program - Controls calls to all routines and main computational loop 


OCEANSUB 


Contains: 

WRITENITRATE,PROFHEADS,WTZPROFS,W_SERIES,W_SOLAR,WRITKT 


LAST TZ SZ 


Outputs temperature and salinity fields at last time step of computation 


LEVMLD 


Computes the mixed-layer depth 


LIDATA 


Opens and reads light data file 


NAVAER 


Computes aerosol parameters according to a simplified version of Navy model 


NH4DOSE 


Computes ammonium dosage for each daily interval 


NH40XTERM 


Computes ammonium to nitrate conversion term 


NH4RIGHT 


Computes source and sink terms at the right hand side of the ammonium equation 


N03RIGHT 


Computes source and sink terms at the right hand side of the nitrate equation 


PLEFTNEWPZ 


Assigns values to tridiagonal matrix that originate from P and Z equations 


PLEFTNEWNN 


Assigns values to tridiagonal matrix that originate from NH 4 and NQ 3 equations 


PRIGHT 


Computes source and sink terms at the right hand side of the phytoplankton equation 


PRIMARY 


Computes gross, net and new production 
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PROFHEADS 


PSTAR 


RD EKMAN DATA 


READPRF 


RETR 


RMODE 


RICH 


RIGHTHANDNN 


RIGHTHANDPZ 


SAVELAST 


SFCRFL 


SHALLOW 


SINKING 


STATE 


STRESS 


SUB 


WATERRAD 


WRITENITRATE 


WRITASC 


WRT DIAG 


WTZPROFS 


ZRIGHT 


Writes header information for the file containing the vertical profiles 


Computes the non-dimensional entrainment 


Contains modules to read, calculate and manipulate Ekman pumping velocit; 


Reads OGCM m, \\ T , and vr fields 


Computes mixed-layer retreat (see RMODE) ___ 


Computes ML retreat while conserving buoyancy, temp., salinity, and potential ener 


Computes all biological quantities and forms right hand side of differential equations 


Computes diffusion and advection terms for the ammonium and nitrate equations 


Computes diffusion and advection terms for the phytoplankton and zooplankton equations 


Save all biological variables at last time step for subsequent hot start of model 


Computes surface reflectance for direct and diffused components separately 


Limits minimum mixed-layer depth to a prescribed value (ex.: 101 cm) 


Computes sinking rate of phytoplankton ____ 


Computes buoyancy from temperature and salinity using simplified UNESCO equation 


Computes sea-surface wind stress ____ 


Contains: 

OPTICS, INITPZN,INIT1JNIT2,INIT3,PRIGHT,ZRIGHT.NH4RIGHT,N03RIGHT, SIN 
KING,PRIMARY,EKMW,COMPUTE_KV 


Computes the solar zenith angle 


This routine is used in combination with ICOND _ 


Computes water column downweiling radiance 


Writes nitrogen production/loss to a file 


Writes time series of bio-physical variables to a file 


Writes each term of the 4-component equations for diagnostic purposes 


Writes daily time series of vertical profiles for selected model variables 


Computes source and sink terms at the righthandside of the zooplankton equation 
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B.2 EC01D with Chen et al. MLM Coupling or OGCM-driven Options 

This section describes the input parameters, nominal values, input and output files, program flowchart, and 
FORTRAN routine names and functions for the ecosystem model using either the Chen et al. MLM or input 
from the OGCM. 

B.2.1 Input Parameters and Nominal Values 

As previously described in section Bl.l, the input parameters and nominal values used by the ecosystem 
model are read in from a file (filename. par) specified in the command line when executing the model run: 

%ecold < filename. par & 

The initial and forcing fields required for the multiyear runs are read in from other files specified in section 
B2.2. The format for all parameters in filename.par is provided in the include file ‘reading.inc’. The 
following table summarizes each of the input parameters contained in the file. 

IE= Upwelling option: 

0: use constant vertical velocity 
1: evoke Ekman upwelling calculation 

CB= Biological model option: 

0: no biological model evoked 
1 : couple mixed layer model with biological model 

IP= Switch for initial nutrient concentration file: 

1 : use default values 

2: read initial nitrate and ammonium concentrations from units 7,8 

3: read initial nitrate and ammonium concentrations from hot start file (7) 

4: read initial nitrate concentration from unit 7 and use default ammonium 
values 

DD= Switch to print results: 

0: do not print results 
1: print results 

DB= Switch to print daily-averaged time series of surface and depth-averaged values: 

0: do not print time series 
1: print time series 

CV= Switch to control writing of hot start file: 

0: do not write hot start file 
1 : write hot start file 

HR = Controls frequency (in hours) of output of model variables (except for profiles) 

HP = Controls frequency (in hours) of output of vertical profiles 

vapr= 2.79 Atmospheric water content (cm) 

uo= 0.34 Ozone concentration (Dobson) 

a_type= Aerosol type: 

0: Continental 
1: Oceanic 


Cldbb=0.50 


Cloud correction coefficient 



Cldaa=0.53 
dt=0.3600D+04 
dz =0.10000+01 
Dep=250.0 
AKD=0.0525 
ZONE=25.0 
BOTKV=8.65 
We= 0.3 
LAT= 0.0 
LON= 165.0 


IW= 


m= 0. 1 
G0= .8511 
Kno3= 1.0000 
Kah4= .5000 
c= 0.02 
Agk= .0633 
Ark= .0633 
Ideep = 400.0 
Ishal= 400.0 
mixed= 60.0 
NH4= 0.1 
N03= 23.0 
P= .1000 
Z= .1000 
pk= 5.0 
Knir= 2.5 
K v = 1.0 
CHLN= 1.0 
DI = 4 
cs= 2 
fc= 1 .0 
Smax= 0.0 
g= 0.04 
Rm= 7.0 
LAMBDA= 0.8 
GAMMA= 0.3 
rz0= 0.019 
krz= 0.08 
az = 0.6 
ap = 0.6 
TH = 0.0 
MOVING= 24 
START= 19830701 
NDAYS=365 


Cloud correction coefficient 
Model time step (seconds) 

Vertical grid resolution (meters) 

Maximum depth for model calculations 

Coefficient for the vertical profile of diffusion for the bio model 
Depth of transition layer (meters) 

Minimum vertical diffusivity Kv (m‘/d) 

Upwelling velocity (m/d) 

Latitude in degrees (North is positive) 

Longitude in degrees (East is positive) 

Control for upwelling velocity: 

1: constant We 

2: read in We from external source 

Phytoplankton death rate, d 1 

Phytoplankton growth rate at 0°C, d' 1 

Half saturation of nitrate uptake, mg-at N m 3 

Half saturation of ammonium uptake, mg-at N m 3 

Ratio of phytoplankton respiration to growth 

Temperature sensitivity of algal growth, 1/°C 

Temperature sensitivity of algal respiration, 1/°C 

Maximum photoacclimation parameter (> 60 m), pmol quanta/nr/s 

Maximum photoacclimation parameter (< 60 m), qmol quanta/nr/s 

Critical depth (meters) for Ishal and Ideep calculation 

Default initial NH4 value (if input file not assigned) 

Default initial N03 value (if input file not assigned) 

Initial Phytoplankton concentration, mg-at N m 3 
Initial Zooplankton concentration, mg-at N m 3 
Nitrate-sensitive coefficient due to ammonium 
Attenuation coefficient for Near Infrared (m ’) 

Minimum value for K v (nr/d) 

Chlorophyll-a:Nitrogen ratio, mg/mg-at N 

Index to control the method to compute surface diffusivity 

Index for choosing the method to compute the diffusivity profile 

Multiplicative factor to apply to vertical diffusivity 

Sinking speed (m/d) 

Zooplankton death rate, d 1 
Zooplankton maximum grazing rate, d 1 
Ivlev constant, m 3 /mg-at N 
Unassimilated zooplankton ingested ratio 
Not used - reformulated 

Temperature sensitivity of zooplankton respiration, 1/°C 
Fraction of dead zooplankton converted to ammonium 
Fraction of dead phytoplankton converted to ammonium 
Minimum C threshold for zooplankton grazing, mgC/m 3 
Moving average window in hours for light energy 
Starting date for calculations, year/month/day 
Total number of days of simulation 


IMLD= 1 


Switch for depth integration limits (Tl) 
1: surface to mixed-layer depth 
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T2= 50.0 
T3= 250.0 
DLEVEL= 250.0 
INTV= 10.0 
Anmax= 02.00 
Dmin= .0095 
KD = 0.0360 
Cpel =0.80 
CN =0.1 
CM =0.85 
GAMA=0.23 
T AUMN=0. 00005 
KVALE=0.0 
HP= 20.0 
Hmin=5.0 
Hmax= 150.0 
WPC=10.0 
EP=0.00001 
KTB=0.0 
KVB=0.0 
RIC1=0.65 
RIC2=0.25 
GAM 1=0.4 
GAM2=0.4 

IHP=0 


COUPLED = 0 
UPWEL=1 


2: surface to bottom 

Integration depth limits 

Integration depth limits 

Maximum depth level for model output 

Depth interval for output 

Maximum rate of nitrification, nmol/day 

Minimum irradiance (300 - 470 nm) for nitrification, W/nr 

Half-saturation dosage for nitrification photoinhibition, W/nT/nm 

Remineralization coefficient for fecal pellets 

Dissipation coefficient of Equation 1.2 

Dissipation coefficient of Equation 1.2 

Coefficient for partial solar radiation (PAR 280 - 700 nm) 

Minimum value of wind stress 
Coefficient of flux adjustment 
Attenuation depth for shortwave radiation 
Minimum mixed layer depth (meters) 

Maximum mixed layer depth (meters) 

wp*wpc 

epsilon 

Background mixing coefficient for temperature and salinity (mVs) 

Background mixing coefficient for u and v (m“/s) 

Parameter for const. Richardson No. profile (instability adjust., J. Price, WHOI) 

Parameter for constant Richardson number profile (after J. Price, WHOI) 

Parameter for constant Richardson number profile (after J. Price, WHOI) 

Parameter for constant Richardson number profile (after J. Price, WHOI) 

Switch for attenuation depth scheme 

0: calculate attenuation depth from light profile 
1: use constant value of 17 meters 

Control index in main module to couple (l)/uncouple(0) bio turbidity in water 
with light absorption coefficient 

Control index to included l)/exclude(0) upwelling in computation 


B.2.2 Input and Output files 

The following logic units (LU) and input/output file names are specified in Tilename.par' after the last 
parameter input described in section B2.1. The file names and data sets that they contain are only given as 
examples; other data sets for different time periods and geographical areas can be specified by the user. 
Note that LU=TMP means temporary logic unit. 

Logic units (LU) for input! output files 

1 2 3 11 12 13 14 7 8 9 80 90 


Input files 


/user_pathname/xbt250.dat 

LU=01; 

/user_pathname/testmet8792.dat 

LU=02; 

/user_pathn ame7raguw8792.dat 

LU=03; 


Initial temperature profiles from XBTs 
Environmental data from TOGA/TAO (1987-1992) 
OGCM vertical velocities (1987-1992) 
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Output Files 


/user_pathname/namef.asc 

LU= 1 1 ; 

/u ser_pathname/ name p.asc 

LU=12; 

/user_pathname/namet.asc 

LU=13; 

/user_pathname/namen.asc 

r 

C 

11 

/user_pathname/diag_bal_name.dat 

LU=88; 


Daily time series of selected diagnostic and predicted 
model variables (see Table A6) 


Daily time series of column-averaged and surface 
values for selected model variables (see Table A8) 

More daily time series of selected diagnostic and 
predicted model variables (see Table A9) 


Diagnostic file containing daily profiles of all equation 
terms (Phyto, Zoo, NCf and NH 4 ) 


Daily time series of vertical profiles for selected model 
variables (see Table A7) 


Input files 

/user_pathname/nitr 1 65e.dat 
/u ser_path n ame/w i nn h4 . d a t 

Output Files 

/user_pathname/lastdyn.bin 

/user_pathname/lastbio.bin 

Input files 

/user_pathname/uvatmodd.dat 

/user_pathname/wv_chl_absp.dat 

/user_pathname/asnh4.dat 

/user_pathname/testrad8792.dat 

/user_pathname/uvtw83792n.dat 

/user_pathname/uvtc83792.day 

/user_pathname/mixed83792.dat 


LU=07; Climatological nitrate (N0 3 ) concentration profile 
LU=08; Climatological ammonium (NH 4 ) concentration profile 

LU=TMP; Binary file containing dynamic variables from last time 
step to be used for hot starting the model 

LU=TMP; Binary file containing biological variables from last time 
step to be used for hot starting the model 

LU-TMP; Spectral absorption coefficients for Gregg-Carder model 

LU=TMP; Attenuation coefficients due to chlorophyll 

LU=TMP; Ammonium nitrification coefficients 

LU=TMP; Radiation data from TOGA-TAO (1987-1992) 

LU=09; Time series of u, v, w and temp from OGCM 

LU=80; Monthly atmospheric u, v\ T , cloud fraction data from 
ISCCP 

LU=90; Mixed layer depth time series (19983-1992) from 
OGCM 


B.2.4 Ecosystem Model Flowchart 

Figure B2.1 shows a flowchart of the second and third versions of the ecosystem model which features 
options for a stand-alone calculation of the temperature and vertical diffusivity profiles using the Chen et al. 
[1994] mixed-layer model, and one-way coupling of velocity and temperature fields derived from the 
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OGCM. The flowchart identifies each main FORTRAN module in a separate box, with the arrows 
indicating input/output data flows and time stepping loops. 
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Figure B2.1. Flowchart of ecosystem model with one-way coupling to OGCM physical parameters. 
Note that there is an option to couple the biology with the Chen MLM. 
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B.2.5 Routine Names and Functions 

To better describe the functions of each FORTRAN module in the flowchart of Figure B2.1, a short 
explanation of each FORTRAN module is provided in Table B2.1. The list of modules in Table B2.1 is 
more extensive than the modules appearing in Figure B2.1 because some FORTRAN modules contain calls 
to subordinate FORTRAN routines and functions. The executable file for the ecosystem model is built by 
compiling and linking all relevant FORTRAN modules using a ‘Makefile’ provided with the source code 
package. 


Table B2.L. Names and functionality of ecosystem model FORTRAN routines. 


Routine Name 

Description 

ABSPAR 

Contains WATERRAD - Computes water column irradiance 

ALGO 

Contains: 

RIGHTHANDPZ,RIGHTHANDNN,PLEFTNEWPZ,PLEFTNEWNN,PLEFT,GAUSS 

ATMODD 

Model of atmospheric transmittance of solar irradiance through marine atmosphere 

BACKZ 

Converts a coordinate to z coordinate 

COMP_ESTAR 

Computes running daily mean of light at each depth (E*) 

CVMIX 

Convective adjustment of T, S, B, U, and V fields for the Chen mixed-layer model 

DAYTIM 

Converts Gregorian to Julian, and Julian to Gregorian dates 

DCHEN 

This module contains routines for the Chen mixed-layer model 

DENS 

Computes seawater density from temperature, salinity and pressure 

DIFFUSION 

Contains mixing and advection modules (BKMIXINGW,BKMIXING,ADVECTION) for 
Chen's MLM 

ENERGY 

Computes energy balance within the mixed layer 

FROUIN 

Computes clear sky solar irradiance using Frouin et al. method 

GAUSS 

Gauss elimination implicit matrix solver for ecosystem model (Crank-Nicholson scheme) 

GREGG 

Computes spectral solar irradiance using the Gregg and Carder method. Contains: 
GR EGG_C ARDER 

GROWTH 

Computes respiration and growth. Contains: 
GROWTH l ,GRO WTH2 ,RESIZ 

HYBRID ID 

Chen's ID MLM using hybrid vertical mixing scheme 

IMPMIX 

Adjusts vertical profiles according to the depth change of the mixed layer 

INIT 

Initializes temperature and salinity fields 


INTERPRF Reads u, i\ T OGCM profiles at 10-m intervals, interpolates values onto a 1-m spacing, 

calculates the Richardson number and mixed-layer depth. Contains: 

RD_UVTW_PROF, INTERPOL, DVARDZ,RICHARDSON,MIXED_LAYER,REALFT 

INITPZN Initializes the phytoplankton, zooplankton, ammonium and nitrate fields 

JPMIX Relieves gradient Richardson number instability using J. Price’s criterion 

KTMIX Calculates vertical mixing using the Kraus-Turner scheme 

KTMODEL Main program - reads/writes data and controls calls to all routines and main computational 

loop 

KTMSUB Contains: 

WRITENITRATE,PROFHEADS,WTZPROFS,W_SERIES,W_SOLAR,WRITKT 

Outputs temperature and salinity fields at last time step of computation 

Opens and reads light data file 

Computes aerosol parameters according to a simplified version of Navy model 

Computes Ammonium dosage for each daily interval 

Computes ammonium to nitrate conversion term 

Computes source and sink terms at the right hand side of the ammonium equation 

Computes source and sink terms at the right hand side of the nitrate equation 

Assigns values to tridiagonal matrix that originate from P and Z equations 

Assigns values to tridiagonal matrix that originate from NH 4 and NQ 3 equations 

Computes source and sink terms at the right hand side of the phytoplankton equation 

Computes gross, net and new production 


LIDATA 

NAVAER 

NH4DOSE 

NH4QXTERM 

NH4RIGHT 

NQ3RIGHT 

PLEFTNEWPZ 

PLEFTNEWNN 

PRIGHT 

PRIMARY 




























































63 


PROFHEADS 

rd,ekman_pata 

REAPPRF 

RICH 

RIGHTHANDNN 

RIGHTHANDPZ 

RICHARDSON 

SAVELAST 

SFCRFL 

STATE 

SUB 


Writes the header for the file containing the time series of vertical profiles 

Contains modules to read, calculate and manipulate Ekman pumping velocity 

Reads OGCM u, v, T , and vr fields 

Computes all biological quantities and forms right hand side of differential equations 

Computes diffusion and advection terms for the ammonium and nitrate equations 

Computes diffusion and advection terms Tor the phytoplankton and zooplankton equations 
Computes the vertical mixing coefficient (K v ) using the Pacanowski and Philander scheme 

Save all biological variables at last time step for subsequent hot start ol model 

Computes surface reflectance for direct and diffused components separately 

Computes buoyancy from temperature and salinity using simplified UNESCO equation 
Contains: 


OPTICS ,INITPZN,INIT1,INIT2,INIT3,PRIGHT,ZRIGHT,NH4RIGHT,N03RIGHT,SINKI 
NG,PRIMARY.EKMW,COMPUTE_KV,DTKV,DTKV8,KD_PROFILE, KVISINSIGMA 


SUNINFO 

Computes the solar zenith angle 

TKEO 

Computes total turbulent kinetic energy to be used in Chen’s MLM 

TRIDIAG 

Tridiauonal matrix solver for the implicit Crank-Nicholson Finite difference method 

WRITENITRATE 

Writes nitrogen production/loss to a file 

.WTZPROFS 

Writes the daily time series of profiles to unit 12 

W_SERIES 

Writes the daily time series of surface and vertically averaged variables to unit 13 

WRT_DIAG 

Writes each term of the 4-component equations for diagnostic purposes 

ZRIGHT 

Computes source and sink terms at the righthandside of the zooplankton equation 
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